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1.0  Studies  of  Neutral  Atmospheric  Densit' 


1 . 1  Introduction 


Three  studies  were  undertaken  during  the  contract  period. 
These  were:  a)  the  Air  Force  Reference  Atmosphere  (AFRA) 
model  evaluation  and  selection  study,  b)  the  Density 
Variability  Study,  c)  revision  of  the  Jacchia  70  Tides,  and 
development  of  the  Air  Force  Reference  Atmosphere  Supplements 
1986.  In  addition,  several  smaller  studies  such  as  the  Total 
Atmospheric  X-ray  Intensity  and  the  Particle/Density  Study 
were  undertaken,  but  are  not  reported  here. 

The  Air  Force  Reference  Atmosphere  model  evaluation  and 
selection  process  began  in  November  1984  by  considering  the 
MS IS  83  and  Jacchia  70  Tides  models.  Numerous  data  bases 
were  prepared  during  the  course  of  this  study,  and  a  packed 
data  base  format  was  developed  and  implemented  in  software 
(See  Section  1.4  below).  By  March  1985,  the  study  had  been 
expanded  to  include  all  the  available  data  bases  at  the  time 
(i.e.  AE/S3-1 ,  S3-4,  SETA-1,  SETA-2,  and  SETA-3).  The  models 
under  consideration  at  this  point  were  the  Jacchia  70, 

Jacchia  70  Tides  (revised),  Jacchia  71,  MSIS  77,  MSIS  79, 

MSIS  83A,  and  MSIS  83B.  The  goal  of  this  effort  was  to 
select  a  basic  model  for  the  reference  atmosphere  in  the  130 
to  200  kilometer  region.  This  work  was  completed  in  April 
1985,  with  the  result  that  two  models,  the  Jacchia  71(1) 
and  the  MSIS  83(2),  were  chosen  for  consideration  in  future 
studies.  The  results  of  this  study  led  into  the  Density 
Variability  Study  and  the  Air  Force  Reference  Atmosphere 
Supplements  1986. 
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The  second  major  work  effort  undertaken  during  the  contract 
period  was  the  Density  Variability  Study.  The  goal  of  the 
Density  Variability  Study  was  to  recommend  an  atmospheric 
density  model  for  operational  use  and  to  provide  error  bounds 
on  the  selected  model.  The  scope  of  the  work  was  outlined  in 
October  1985,  and  the  main  objectives  were: 

i)  Quantify  accelerometer  measurement  accuracy  vs.  altitude 

ii)  Estimate  accuracy  of  selected  density  models 

iii )  Study  spatial  correlation  of  density  errors  during 
quiet  and  disturbed  geomagnetic  conditions. 

A  schedule  was  agreed  on  and  the  bulk  of  the  work  was 
completed  in  December  1985,  continuing  into  early  1986.  A 
first  draft<3)  of  the  Density  Variability  Study  was  reviewed 
and  commented  upon  in  June  1986. 

Two  considerations  are  worth  noting.  First,  it  was  suspected 
that,  due  to  the  manner  in  which  the  data  were  collected,  the 
resulting  estimates  of  atmospheric  density  are  correlated  in 
time  (auto-correlated).  A  study  was  undertaken  (See  Section 
1.6)  to  investigate  this  hypothesis,  and  this  conclusion  was 
confirmed.  While  this  effect  was  shown  not  to  affect  the 
standard  deviations  as  currently  calculated  by  the  STAT 
program,  it  would  affect  the  calculation  of  any  error  bounds 
on  predicted  densities.  Secondly,  in  January  1986,  it  was 
noted  that  there  was  an  unresolved  question  as  to  the  high 
variability  in  the  measured  to  model  ratios  in  the  AE  data, 
as  indicated  by  the  standard  deviations,  compared  to  the 
variability  obtained  for  the  SETA  data  (See  Table  1.1). 

Hence,  a  study  was  made  using  the  AE-E  data  base  to  see  if 
elimination  of  the  local  time  variability  would  bring  the 
overall  variances  into  agreement.  The  results  of  this  study 
are  included  here  in  Section  1.7.  This  analysis  shows  that 
although  there  are  significant  unmodeled  local  time 
variations  in  the  AE-E  data  base  for  both  the  MSIS  83  and 
Jacchia  71  models  in  a  statistical  sense,  the  resulting 
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overall  standard  deviation  still  remains  at  around  12%,  well 
above  the  errors  for  the  SETA  satellites.  Further  investi¬ 
gative  work  along  the  lines  of  the  analysis  in  Section  1.6 
might  help  to  resolve  this  question. 

The  most  recent  effort  undertaken  was  the  Air  Force  Reference 
Atmosphere  Supplements  1986.  Work  began  on  this  effort  in 
January  1986  when  the  files  containing  the  Air  Force 
Reference  Atmosphere  in  the  70-120  kilometer  region  were 
received(4) .  Preliminary  comparisons  of  the  "Forbes"  model 
with  the  MSIS  83  model  were  made.  An  attempt  was  made  to 
develop  a  model  which  interpolated  between  the  Forbes  and  the 
MSIS  83  models  in  the  104-120  km  region  while  simultaneously 
satisfying  hydrostatic  equilibrium.  When  this  attempt  failed 
to  produce  acceptable  results,  and  after  consultation  with 
AFGL  personnel,  the  decision  was  made  to  drop  the  hydrostatic 
equilibrium  requirement  and  to  use  cubic  spline  fits  to 
interpolate  in  the  104-120  km  region  between  the  "Forbes" 
model  and  MSIS  83.  Work  on  the  Air  Force  Reference 
Supplements  1986  has  since  been  completed  and  the  results 
have  been  delivered  to  the  initiator.  These  results 
consisted  of  110  pen  plots  and  756  pages  of  tables. 


Table  1.1  Statistical  Properties  of  Total  Mass  Density 


Loa  Ratios 

to  MS  IS 

83  Model; 

:  Data  to 

200  km 

Data  Set 

N 

X 

s 

^2. 

AE-C 

62,044 

0.094 

0.152 

-0.619 

8.156 

AE-D 

29,762 

-0.029 

0.157 

-0.775 

10.857 

AE-E 

35,629 

0.014 

0.132 

-0.457 

8.812 

S3-1 

27,315 

0.050 

0.156 

-0.225 

12.311 

SETA-1 

20,923 

-0.073 

0.087 

-0.718 

9.946 

SETA-2 

105,809 

-0.150 

0.107 

-0.503 

5.353 

SETA-3 

( Jul-Nov) 

652,238 

-0.144 

0.092 

0.192 

5.508 

N  is  the 

sample  size 

i,  and  x. 

s,  Jbw 

and  b2  (defined 

below ) 

are  the  first  four  sample  moments  about  the  sample  mean. 


Table  1.2  Neutral  Atmospheric  Density  Satellite 
Experiment  Flight  Histories 


SATELLITE  JULIAN 

SEASON 

ALTITUDE 

INCL 

MEAN 

# 

# 

# 

SECS 

DATE 

RANGE  (KM) 

F10.7 

FILES  RECS 

POINTS 

CASE 

AE-C 

73353 

74352 

FWSS 

135 

245 

68.4 

86.9 

62044 

44 

S3-1 

74309 

75120 

FWS 

135 

245 

97 

76.3 

27315 

44 

AE-D 

75281 

76029 

FW 

140 

245 

90 

75.8 

29762 

44 

AE-E 

75325 

76322 

FWSS 

135 

245 

20 

73.3 

35629 

44 

AE/S3-1 

73353 

76322 

135 

245 

23 

12929 

154750 

44 

S3-4 

78089 

78223 

SS 

162 

285 

96* 

145.8 

161 

4756 

54081 

26 

SETA-1 

79079 

79100 

SPRING 

168 

249 

96* 

189.8 

1 

5222 

62649 

23 

SETA-2 

82137 

82332 

SSF 

168 

297 

96* 

170.0 

114 

37148 

443805 

110 

SETA-3 

83201 

84075 

SFWS 

163 

385 

96* 

116.0 

150 

246624 

2956869 

655 

S85-1 

84177 

84283 

SF 

179 

257 

96* 

81.3 

51 

20685 

247330 

113 

*  SUN  SYNCHRONOUS  ORBIT 


1.2.1  Overiew 


The  STAT  program  is  the  primary  analysis  tool  used  in 
studying  the  neutral  atmospheric  density  data  base.  STAT 
typically  breaks  the  data  down  into  a  set  of  bins,  computes, 
and  prints  a  report  consisting  of  the  means  and  standard 
deviations  in  each  bin  for  selected  satellite  and  model 
combinations.  The  data  base  consists  of  estimates  of 
atmospheric  density  computed  from  accelerometer  measurements 
of  total  satellite  drag.<5'6)  A  list  of  the  satellites  and 
their  flight  histories  which  make  up  the  data  base  is  given 
in  Table  1.2.  The  current  list  of  data  base  tapes  in  machine 
readable  form  is  given  in  Table  1.3.  The  data  base  format  is 
given  in  Table  1.4.  The  STAT  program,  based  on  the  input, 
requests  the  appropriate  data  base  tape  (or  tapes)  via  CALLs 
to  PACKLIB  routines  (See  below.),  and  calculates  and  prints 
the  requested  statistics  for  the  requested 
satellite( s )/model( s )  combination.  Currently,  the  only 
restriction  is  that  all  the  requested  models  (up  to  a  maximim 
of  four)  must  be  stored  on  the  same  data  base  tape.  A 
flowchart  of  the  STAT  program  is  given  in  Figure  1.1. 


Table  1.3  Air  Force  Reference  Atmosphere 
Neutral  Atmospheric  Density  Data  Base 
(File;  DBINDEX/UN»BRYANTC) 


f 

f 

t 

t 

} 

4 

i 

I 

3 


l 


CC4494 

H  SETA-1 

22371 

23370T 

24MSIS86 

26MSIS83B 

i 

CC1465GE 

H  AE/S3-1 

22J71 

23370T 

24MSIS77 

26MSIS83B 

■ 

CC5180GE 

H  AE/S3-1 

22J71 

23370T 

24MSIS77 

26MSIS83B 

CC0594GE 

H  AE/S3-1 

22371 

2337 OT 

24MSIS77 

26MSIS83B 

CC0394 

H  S3-4 

22J71 

23370T 

24MSIS77 

26MSIS83B 

» 

CC5061 

H  SETA-1 

22J71 

23370T 

24MSIS77 

26MSIS83B 

> 

CC0764 

PH  SETA-2 

22J71 

23370T 

24MSIS77 

26MSIS83B 

■ 

'i 

CC2462GE 

PH  SETA-3 

22371 

23370T 

24MSIS77 

26MSIS83B 

i 

CC1403GE 

PH  SETA-3 

22371 

23370T 

24MSIS77 

26MSIS83B 

% 

m 

CC4274 

PH  S85-1 

22371 

23370T 

24MSIS77 

26MSIS83B 

CC2478 

PH  AE/S3-1 

24MSIS77 

26MSIS78/79 

V 

( 

0S0475  ; 

SI  S3-4 

23377 

24MSIS77 

26MSIS78/79 

i 

CC4953  : 

SI  H  SETA-1 

23377 

24MSIS77 

26MSIS78/79 

i 

% 

CC0221 

PH  SETA-2 

23377 

24MSIS77 

26MSIS78/79 

s 

CC2466GE 

H  AE/S3-1 

23377 

24MSIS83A 

26MSIS83B 

4 

■ 

OS0637 

SI  H  S3-4 

24MSIS83A 

26MSIS83B 

| 

CC4954 

SI  H  SETA-1 

24MSIS83A 

26MSIS83B 

5 

CC2125 

PH  SETA-2 

24MSIS83A 

26MSIS83B 

• 

V 

CC2484GE 

H  AE/S3-1 

24370T 

\ 

CC0413 

H  S3-4 

24370T 

j 

CC1065 

H  SETA-1 

24370T 

CC2165 

PH  SETA-2 

24370T 

5 

CC2493GE 

<AE/S3-1 

22364 

23US66 

24371 

38USSR 

! 

OS0472 

SI  S3-4 

22364 

23US66 

24371 

38USSR 

« 

OS0492 

H  SETA-1 

22364 

23US66 

24371 

38USSR 

0 

CC0061 

PH  SETA- 2 

22364 

23US66 

24371 

38USSR 

a 

a 

CC2515GE 

<AE/S3-1 

22JWB 

23L/N 

38US62 

a 

OS0473 

H  S3-4 

223WB 

23L/N 

24US62 

* 

0S0494 

H  SETA-1 

22JWB 

23L/N 

24US62 

CC0081 

PH  SETA- 2 

223WB 

23L/N 

24US62 

CC2523GE 

<AE/S3-1 

23DENSEL 

24370 

38373 

i 

OS0474 

SI  S3-4 

22DENSEL 

23370 

24373 

) 

0S0496 

H  SETA-1 

22DENSEL 

23370 

24373 

i 

CC0101 

PH  SETA-2 

22DENSEL 

23370 

24373 

• 

CC5019 

SETA-1 

32371KP-1 

OS0098 

SI  SETA-1 

32371KP-1 

1 

CC0583 

P  SETA- 2 

32371KP-1 

« 

CC5108 

P  S85-1 

32371KP-1 

3 

Table  1.4  Neutral  Atmospheric  Density  Data  Base  Format 


HEADER  RECORD 

WORD#  CONTENTS 

0.1  WORD  COUNT  (40) 

0.2  GROUP  COUNT  (1) 

1  SATELLITE  ID 

2  EXPERIMENT  NAME 

3  ALTITUDE  (KM) 

4  100*( NUMBER  OF  FILES)  +  (NUMBER  FOR  THIS  FILE) 

5  DATA  BASE  CREATION  DATE  (YYDDD) 

6-40  BLANK 


DATA  RECORD 
WORD# 

~on 

0.2 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 


CONTENTS 
WORD  COUNT  (40) 

GROUP  COUNT  (12) 

ORBIT  NUMBER 

JULIAN  DATE  OF  DATA  (YYDDD) 

GREENWICH  TIME  (SECONDS) 

GREENWICH  HOURS 

GREENWICH  MINUTES 

GREENWICH  SECONDS 

LOCAL  TIME  HOURS 

LOCAL  TIME  MINUTES 

LOCAL  TIME  SECONDS 

UPLEG/DOWNLEG 

DAY/NIGHT 

SPUN/DESPUN 

GEOGRAPHIC  LATITUDE 

GEOGRAPHIC  LONGITUDE  (EAST) 

GEOMAGNETIC  LATITUDE 

GEOMAGNETIC  LONGITUDE  (EAST) 

GEOMAGNETIC  LOCAL  TIME 

MEASURED  DENSITY  (G/CM**3) 

J71  MODEL  DENSITY 

J71  MODEL  DENSITY  (KP*2) 

MSIS  MODEL  DENSITY 
MEASURED/ J71  RATIO 
MEASURED/ J71 ( KP=2J  RATIO 
MEASURED/MSIS  83  RATIO 
AP  (DAILY  AVERAGE) 

AP  (6.7  HOUR  LAG) 

KP  (6.7  HOUR  LAG) 

F10.7  1  DAY  LAG] 

F10.7  (81  DAY  AVERAGE) 
ALTITUDE  (KM) 

KP  (3  HOUR  LAG) 

KP  (6  HOUR  LAG) 

KP  (9  HOUR  LAG] 

KP  (12  HOUR  LAG) 

AVERAGE  OF  WORDS  31-34 

J77  MODEL  DENSITY 
MEASURED? J7 7  RATIO 
EMPIRICAL  MODEL  RATIO 
MEASURED/EMPIRICAL  RATIO 


SBINS  VARIABLE 
- BTN - 


'ORBIT' 
'  JDATE ’ 
'UT' 


'LT' 


'SD' 

’  GGLAT ' 

’  GGLONG ' 
'  GMLAT ' 

'  GMLONG ' 
'  GMLT ’ 


'AP' 

' AP6 . 7 ' 
'KP' 

' F10. 7 ' 

'F10.7BAR' 

'ALT' 

'  KP3 ' 

'  KP6 ' 

'  KP9 ' 

' KP12 ' 

' KPAVG ’ 
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Figure  1 . 1  STAT  Flowchart 


1.2.2  STAT  Input  Description 


The  STAT  program  is  executed  via  the  STAT  procedure  which 
is  stored  under  UN=BRYANTC.  To  get  a  copy  of  the  STAT 
procedure,  say  "GET, STAT/UN=BRYANTC"  at  the  command  level 
while  logged  in  under  NOS  2.  Parameters  for  the  STAT 


procedure  are: 

TABLES 

File  name  under  which  to  store  the  printed 
tables.  (May  be  omitted. ) 

LASER 

YES  or  NO.  Send  tables  to  the  the  Xerox-2700 
laser  printer  or  not.  Default  is  NO. 

COPIES 

Number  of  copies.  Default  is  one  copy. 

TITLE 

Title  for  laser  output  (1-10  characters). 

SUATEK 

File  name  under  which  to  store  SUATEK  format 
file  (direct  access).  Default  is  SUATEK. 

BINS 

YES  or  NO.  Print  bin  structure  or  not. 
Default  is  NO. 

CPUSAGE 

YES  or  NO.  Print  cp  usage  statistics  or  not. 
Default  is  YES. 

PACKED 

YES,  NO,  or  UNDEF.  Packed  data  base  or  not, 
if  known.  Default  is  UNDEF. 

UN 

User  name  for  laser  output  banner.  No 
default . 

Input  for  a  case  consists  of  a  single  title  card,  followed  by 
the  NAMELIST  SINPUT  input,  followed  by  the  NAMELIST  SBINS 
input.  The  user  may  supply  a  subroutine  named  DATAMOD  as  a 
second  record  of  input.  If  the  user  sets  the  variable 
USERSUB  true  in  conjunction  with  this,  then  the  subroutine 
DATAMOD  will  be  CALLed  on  every  pass  through  the  main 
processing  loop  after  checks  on  altitude,  date  range  and 
spun/despun  have  been  made.  This  gives  the  user  freedom  to 
test  and/or  modify  each  data  point  prior  to  further 
processing.  A  sample  DATAMOD  subroutine  is  given  in  Figure 
1.2.  Cases  may  be  stacked  to  any  degree  desired,  limited 
only  by  the  amount  of  cp  time  requested  on  the  JOB  card. 


IF  (J  .EQ.  0)  THEN 

CALL  TRANS  ( IHEADER, 40, 1 , CHEADER ) 


CODE  TO  PROCESS  HEADER  OR  SKIP  TO  NEXT  HEADER  RECORD 
SKIP  =  .FALSE. 


ELSE 

CALL  TRANS  ( IDATA( 1 , J ) , 40 , 1 , CDATA( 1 ) ) 


CODE  TO  PROCESS  OR  SKIP  J ’ TH  DATA  POINT 
SKIP  =  .FALSE. 

END  IF 
RETURN 


END 
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•  Title  card  -  80  columns  of  case  descriptive  character  data. 

•  NAMELIST  SINPUT  - 


The  NAMELIST  SINPUT  describes  all  of  the  relevant  information 
for  a  case  except  for  the  description  of  the  bin  structure 
(and  the  title  card). 


SINPUT 


Variable 

Size 

Type 

Description 

Default 

SAT 

9 

CHAR*10 

satellite 

None 

( reguired ) 
MODEL 

4 

INTEGER 

model  number  (1-4) 

None 

WORD 

4 

INTEGER 

word  number 

None 

MODNAME 

4 

CHAR* 10 

model  name 

Predefined 

MODABBR 

4 

CHAR* 10 

model  abbreviation 

None 

The  user  should  specify  either  MODEL,  WORD,  or  MODDABBR. 
MODABBR  is  usually  the  most  convenient  when  doing  statistics 
on  model  ratios.  (See  Table  1.3  for  model  abbreviations.) 
WORD  gives  the  flexibility  to  do  statistics  on  any  variable 
in  the  data  oase.  (See  Table  1.4) 


Variable 

Size 

Type 

Description 

Default 

MODHEAD 
( optional ) 

mo5fmt 

( optional ) 

4 

4 

CHAR* 16 

CHAR*20 

heading  for  model  column 

format  under  which  to 
print  model  statistics 

Predefined 

Predefined 

MRSTATS 

1 

LOGICAL 

true  if  model  ratio 
statistics  are  being 
calculated,  false  o.w. 

ALLDATA 

1 

LOGICAL 

true  if  "all  data"  case 
is  to  be  calculated 

TRUE 

ALT 

2 

REAL 

altitude  range  min  &  max 

0-600 

DAYNITE 

1 

LOGICAL 

true  for  day/night 
partioning,  false  o.w. 

FALSE 

GMLT 

1 

LOGICAL 

true  for  geomagnetic 
local  time,  if  false, 
solar  local  time  is  used 

DAYLT 

2 

INTEGER 

local  time  range  for 
day  ( hhmm  format ) 

0800- 

1800 

NITELT 

2 

INTEGER 

local  time  range  for 
night  ( hhmm  format ) 

1800- 

2400 

ILT 

1 

LOGICAL 

true  if  local  time  is  to 
be  calculated  (packed 
data  bases  only) 

IGMCO 

1 

LOGICAL 

true  if  geomagnetic 
coordinates  are  to  be 
calculated  (packed  data 
bases  only) 

I  GMLT 

1 

LOGICAL 

true  if  geomagnetic  local 
lime  is  to  be  calculated 
(packed  data  bases  only) 
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Variable 

Size 

Type 

Description 

Default 

YYDDD 

2 

INTEGER 

date  range  (yyddd  format) 

None 

YRMODA 

2 

INTEGER 

date  range 
(yymmdd  format) 

None 

DATES 

50 

INTEGER 

set  of  up  to  50  dates  in 
yyddd  format 

None 

SPUN 

1 

LOGICAL 

true  if  spun  data  only 
(AE  satellites  only) 

FALSE 

DESPUN 

1 

LOGICAL 

true  if  despun  data  only 
(AE  satellites  only) 

FALSE 

MEANKP 

1 

LOGICAL 

calculate  and  print 
average  Kp  in  each  bin 

TRUE 

MEANF 

1 

LOGICAL 

calculate  and  print 
average  F10  7  in  each 
bin 

TRUE 

MEANZ 

1 

LOGICAL 

calculate  and  print 
average  altitude  in  each 
bin 

TRUE 

MEANLT 

1 

LOGICAL 

calculate  and  print 
average  local  time  in 
each  bin 

TRUE 

MEANUT 

1 

LOGICAL 

calculate  and  print 
average  universal  time 
in  each  bin 

TRUE 

MEANALL 

1 

LOGICAL 

calculate  and  print 
average  of  each  of  the  5 
variables,  Kp,  F10  7, 
altitude,  local  rime,  and 
universal  time  in  each 
bin 

TRUE 

MEAN 

1 

LOGICAL 

calculate  and  print 
averages  of  the  dependent 
variables 

TRUE 

SD 

1 

LOGICAL 

calculate  and  print 
standard  deviations  of 
the  dependent  variables 

TRUE 

MEANCL 

1 

LOGICAL 

calculate  and  print 
confidence  limits  on  the 
mean  of  the  dependent 
variables 

FALSE 

SDCL 

1 

LOGICAL 

calculate  and  print 
confidence  limits  on  the 
standard  deviation  of  the 
dependent  variables 

TRUE 

ANOVA 

1 

LOGICAL 

Do  a  one-way  analysis  of 
variance  on  the  dependent 
variable 

FALSE 

ALPHA 

1 

LOGICAL 

significance  level  for 
F-rest  in  one-way  ANOVA 

.05 

Default 


Variable 

Size 

Type 

Description 

Defauli 

LNRATIO 

1 

LOGICAL 

calculate  statistics 
using  natural  logarithms 
of  dependent  variables 

TRUE 

LINEAR 

1 

LOGICAL 

calculate  statistics 
using  dependent  variable 
directly 

FALSE 

PCTSD 

1 

LOGICAL 

print  standard  deviation 
as  a  percentage  of  the 
mean 

FALSE 

SUATEK 

1 

LOGICAL 

Write  a  SUATEK  format 
file  and  save  it  as  a 
direct  access  file 

FALSE 

SUABINS 

1 

CHAR* 10 

'MINMAX':  write  bin  mins 
&  maxs  to  SUATEK  file; 
'AVERAGE':  write  bin 
averages  to  SUATEK  file. 

MINMAX 

SUAMEAN 

1 

LOGICAL 

Write  the  mean  Kd, 

F,0  7,  altitude,  local 
time,  and  universal  time 
to  the  SUATEK  format  file 

FALSE 

SUACL 

1 

LOGICAL 

Write  confidence  limits 
on  the  mean  to  the  SUATEK 
format  file 

FALSE 

USERSUB 

1 

LOGICAL 

CALL  the  user  supplied 
subroutine  DATAMDd  in  the 
main  processing  loop 

FALSE 

RAWDATA 

1 

LOGICAL 

true  for  "raw  data" 
print  out,  false  o.w. 

FALSE 

NLPAGE 

1 

INTEGER 

number  of  lines  per  page 
(18<NLPAGE<48,  not 
including  headings ) 

36 

RESET 

1 

LOGICAL 

Reset  the  $ INPUT 
variables  to  their 
defaults  at  the  beginning 
of  the  next  case 

FALSE 

Character  input  must  be  enclosed  in  single  quotes  (').  The 
pairs  of  variables,  LNRATIO  and  LINEAR,  and  SD  and  PCTSD  are 
logical  complements  of  each  other. 


•  NAMELIST  SBINS 


The  user  may  specify  up  to  three  (3)  separate  binning 
variables.  These,  along  with  day/night  partitioning  allow 
for  a  up  to  four-way  breakdown  of  the  data.  (In  the 
following,  "I"  refers  to  the  bin  number,  1,  2,  or  3 . ) 


SBINS 

Variable  Size  Type  Description  Default 


The  user  may  specify  a  particular  bin  variable  using  one  of 
the  following  two  variables. 


BIN  3  CHAR*10  abbreviated  bin  variable 

name  (See  Table  1.4) 

BINVAR  3  INTEGER  bin  variable  number 

(See  Table  1.4) 

The  user  has  a  choice  of  one  of  the  following  six  methods  for 
specifying  each  bin.  (KP2,  KP3,  and  KP4  count  a  one  method. ) 


FROM 

3 

REAL 

These  three  variables 

TO 

3 

REAL 

together  specify  the 

BY 

3 

REAL 

the  bin  structure 
for  the  bin  variable 
named  by  BIN( I )  OR 

BINVAR( I ) 

BMINMAX 

2,36,3 

REAL 

Explicit  enumeration  of 
bin  structure  for  bin  I 

KP2 

1 

INTEGER 

KP2=I  selects  standard 

2  Kp  bins,  (0  to  3+ 
and  4-  to  9)  for  bin 
variable  I 

KP3 

1 

INTEGER 

KP3=I  selects  standard 

3  Kp  bins  (0  to  3,  3+  to 
4+,  and  5-  to  9)  for  bin 
variable  I 

KP4 

1 

INTEGER 

KP4=I  selects  standard 

4  Kp  bins  (0  to  1  +  ,  2-  to 
3+ ,  4-  tO  5+ ,  and  6-  to 

9)  for  bin  variable  I 

MONTH 

1 

INTEGER 

M0NTH=I  selects  monthly 
binning  for  bin  I 

SEASON 

1 

INTEGER 

SEAS0N= I  selects  seasonal 
binning  for  bin  I 

ORBIT 

2000 

INTEGER 

set  of  up  to  2000  orbit 
numbers 

(Only  one  bin  is  allowed  when  ORBIT  is  specified. ) 


Variable 

Size 

TYPe 

Description  Default 

BINNAME 
(optional ) 

40 

CHAR* 14 

bin  name  for  column  Predefined 

heading 

BINFMT 
(optional ) 

40 

CHAR *20 

format  under  which  to  Predefined 

print  out  bin  variable 

1LT 

(optional ) 

1 

LOGICAL 

true  if  local  time  is  to 
be  calculated  (packed 
data  bases  only) 

IGMCO 
(optional ) 

1 

LOGICAL 

true  if  geomagnetic 
coordinates  are  to  be 
calculated  (packed  data 
bases  only) 

IGMLT 
( optional ) 

1 

LOGICAL 

true  if  geomagnetic  local 
time  is  to  be  calculated 
(packed  data  bases  only) 

RESET 
(optional ) 

1 

LOGICAL 

Reset  the  $BINS 
variables  to  their 
defaults  at  the  beginning 
of  the  next  case 

1.3.1  Data  Base  Construction 

The  construction  of  the  neutral  atmospheric  density  data 
bases  described  in  Tables  1.3  and  1.4  is  accomplished  with 
updated  versions  of  programs  FOURMOD  and  DENDB.(7)  Both 
these  programs  have  been  updated  to  produce  on  option  the 
packed  versions  of  the  data  bases  shown  in  those  tables. 

This  format  is  produced  and  accessed  with  the  PACKLIB  library 
to  be  discussed  in  Section  1.4.  In  summary  it  consists  of  15- 
bit  integer  data  words  packed  4  per  full  CDC  60-bit  word. 

For  completeness  we  provide  the  detailed  contents  of  these 
data  bases  in  Table  1.5.  Note  that  only  a  subset  of  the  full 
unpacked  data  base  is  contained.  The  remaining  quantities 
are  easily  derived.  The  implementation  of  this  packed  data 
base  structure  was  motivated  by  the  increasing  amounts  of 
data  to  be  handled.  Program  DENDB  has  also  been  modified 
from  time  to  time  by  substitution  of  density  computation 
modules  for  computation  of  models  other  than  those 
present  in  the  original  package. 

The  complexity  of  the  thermospheric  response  to  geomagnetic 
activity  has  motivated  special  studies  which  have  required 
special  processing.  Program  MRDB  was  written  to  produce 
data  bases  containing  the  ratios  of  the  measured  density  to 
the  Jacchia  71  model  at  fixed  Kp=l.  These  ratios  contain 
the  unattennuated  density  variations  related  to  geomagnetic 
activity,  while  the  other  variations  (altitude,  latitude, 
solar  flux,  etc. )  are  "damped  out"  by  the  approximately 
correct  variations  contained  in  the  Jacchia  71  model. 


MRDB  produces  both  packed  and  unpacked  data  base  formats 
similar  to  those  described  in  Tables  1.3-1. 5,  with  the 
addition  of  the  Kp  index  at  lags  of  -3  hr  to  +15  hr  to  permit 
study  of  density  variations  versus  lag  time  relative  to  Kp 
variations.  Program  OADB  uses  these  data  bases  to  produced 
averaged  ratio  per  latitude  bin  per  orbit,  and  program 
PERCENT  computes  the  percent  changes  of  these  averaged  ratios 
during  period  of  high  geomagnetic  activity  relative  to  the 
values  found  during  the  preceding  quiet  period.  All  these 
results  are  output  onto  data  files  suitable  for  display  by 
program  SUATEK.(8)  Display  of  these  results,  along  with 
proposed  new  geomagnetic  activity  indices,  such  as  the  Total 
Atmospheric  X-ray  Intensity  (TAXI)  index,  permitted 
evaluation  of  these  indices  as  possible  improvements  over  Kp 
for  use  in  density  modelling. 


Header  Record  1:  Same  as  for  unpacked  data  bases. 

Header  Record  2:  Unpacked  offsets  and  scale  values  to  be 

used  to  convert  the  packed  data  words  to  the 
equivalent  60  bit  real  unpacked  values. 

Data  records:  15-bit  data  words,  packed  4  per  60-bit  full 

word;  these  are  preceded  by  a  real  60  bit  word 
containing  the  time  in  the  form  YYDDD.XXXXX 
where  YY  =  the  last  two  digits  of  the  year; 

DDD  =  the  day  number  in  the  year; 

XXXXX  =  fraction  of  the  day  from  UT 
midnight . 

The  15-bit  packed  data  words  are: 

Variable 
Orbit  Number 
Latitude  (°  N) 

Longitude  (°  E) 

Altitude  (km) 

^10.7 

Fi0  7  smoothed 

Kp 

Measured  density  (g/cc) 
Measured/model  density  ratios 


Word  # 
1 
2 

3 

4 

5 

6 

7 

8 

9-12 


1.3.2  PULL/CONSOL/MERGE 


Three  software  programs,  PULL,  CONSOL,  and  MERGE,  have  been 
developed  to  consolidate  selected  models,  originally  residing 
on  different  tapes  onto  a  single  tape.  PULL  "pulls"  the 
data-to-model  ratios  for  selected  models  from  a  full  data 
base  and  packs  them,  time-tagged,  onto  a  temporary  disk  file. 
CONSOL  consolidates  the  model  ratios  on  two  or  three  files 
into  one.  MERGE  merges  the  ratios  residing  on  a  packed 
temporary  disk  into  a  full  data  base.  A  flowchart  of  the 
PULL/CONSOL/MERGE  operation  is  given  in  Figure  1.3. 


1.3.3  HISTO 


HISTO  plots  histograms  of  the  number  of  data  points  versus 
Kp,  solar  flux,  altitude,  geographic  latitude,  and  solar 
local  time  for  a  given  satellite.  Sample  histograms  are 
shown  in  Figures  1.4-1. 8.  Such  histograms  are  useful  in 
showing  the  distribution  of  the  data. 


1.3.4  FREQ /ED I ST 

FREQ  plots  the  relative  frequency  distribution  as  a  function 
of  model  or  ln( model  ratio)  and  superimposes  the  normal 
frequency  distribution  based  on  the  sample  mean  and  standard 
deviation.  The  relative  frequency,  fA,  for  a  given  interval 
is  defined  as  the  fraction  of  the  total  sample  falling  within 
that  interval.  The  first  four  sample  moments,  the  mean, 
standard  deviation,  skewness,  and  kurtosis  are  calculated 
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Figure  1.8  Histogram  of  SETA-1  vs  Local  Time 
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and  displayed  in  a  corner  of  the  plot.  Sample  plots  from  the 
FREQ  program  are  shown  in  Figures  1.9  and  1.10. 


The  mean,  standard  deviation,  skewness,  and  kurtosis  are 
defined  as  follows  <9): 
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n 

iV x* 

-  jO/s) 

kurtosis 

= 

(1/n) 

i?V(  Xi 

-  *T7s)‘ 

ED1ST  plots  the  empirical  distribution  function,  Sn(x),  for 
selected  satellite/density  model  combinations.  The  empirical 
distribution  function,  Sn(x),  is  defined  as  (10> : 

0.  if  x  <  X(1>, 

Sn(  x  )  =  k/n,  if  X(k)  —  x  k=l ,  .  .  .  ,  n-1 , 

l,  it  x  _  X(n) 

where  Xtl) ,  .  .  . ,  X(n)  denotes  the  order 
statistics  of  a  sample  of  size  n. 

The  empirical  distribution  function  is  plotted  on  an  inverse 
normal  probability  scale  versus  a  standardized  abscissa  so 
that  normally  distributed  data  plot  as  a  straight  line.  Let 


(X) 

=  1/42v<j  /  exp[-(t 

-« 

-  m  )/<r  ]2dt. 

let  u 

=  x. 

and  <r 

=  s . 

Then  'l(x)  is  plotted  against  z  =  (x 

-  x)/s,  a  standardized 

deviate.  The  x  and  y  axes  are  labeled  in  terms  of  the 
original  scale  of  the  variable  of  interest  and  cumulative 
probability,  respectively . (UI  Sample  EDIST  plots  are  shown 
in  Figures  1.11  and  1.12. 
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Figure  1.12  SETA-1  Empirical  Distribution  Function  of 
Percent  Difference  from  Jacchia  71  Ratio 
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1.3.5  KLL 


KLL  generates  Kp/latitude/local  time  plots  from  SUATEK  format 
files  produced  by  the  STAT  program.  The  interesting  feature 
of  the  KLL  plots  is  that  the  abscissa  is  non-monotonic.  The 
x-axis,  from  left  to  right,  spans  the  dayside  latitude  range 
-90  to  90  degrees  and  continues  over  into  the  nightside  +90 
to  -90  degrees  latitude  with  appropriate  axis  labeling.  The 
four  Kp  ranges  are  plotted  with  different  symbols.  Sample 
KLL  plots  are  shown  in  Figures  1.13  and  1.14. 


1.4  PACKLIB 


PACKLIB  is  a  collection  of  Radex  written  routines,  originally 
developed  to  read  the  "packed"  data  base  format  tapes. 

PACKLIB  has  since  been  enhanced  and  added  to,  so  that  the 
relevant  routines  are  capable  of  reading  both  the  packed  and 
unpacked  data  base  format  using  a  handful  of  user  friendly 
CALLs.  The  standard  CALLing  sequence  for  processing 
atmospheric  density  data  base  tapes  is  shown  in  Figure  1.15. 

A  number  of  other  routines  useful  in  routine  programming 
applications  have  been  included  in  PACKLIB  as  a  convenient 
repository.  A  list  of  the  PACKLIB  routines,  along  with  their 
argument  lists  is  shown  in  Figure  1.16. 
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SETA-1  SATELLITE  ACCELEROMETER  DATA 
MEAN  LN  (RATIO)*  JACCHIA  71  MODEL 
•  ALL  DATA.  MAX.  ALT.  =  200  KM. 


SETA-1  SATELLITE  ACCELEROMETER  DATA 


SAMPLE  PROGRAM  ILLUSTRATING  THE  USUAL  SEQUENCE  OF  CALLS  TO 
PACKLIB  ROUTINES  WHEN  PROCESSING  ATMOSPHERIC  DENSITY  DATA  BASES. 


INTEGER  DBUNIT 

LOGICAL  SELECT, ILT, IGMCO, IGMLT 

CHARACTER* 10  SAT, MODABBR , MODNAME, CHEADER, CDATA 

COMMON  /HEADER  /  HEADER (40) 

COMMON  /DATA  /  DATA(40,12) 

DIMENSION  MODABBR ( 4 ) , MRWORD ( 4 ) , MODNAME ( 4 ) 

DIMENSION  IHEADER( 40) , CHEADER ( 40 ) 

DIMENSION  I DATA ( 40,12), CDATA ( 40,12) 

EQUIVALENCE  ( HEADER( 1 ) , IHEADER( 1 ) ) , ( DATA( 1,1), IDATA( 1 , 1 ) ) 


DATA  ILT/ . TRUE . / ,  IGMCO/ . FALSE ./ ,  IGMLT/ . FALSE . / 


DBUNIT  =  1 

SAT  =  'SETA-1' 

NMODELS  =  2 
MODABBR ( 1 )  =  ' MSIS83 ' 

MODABBR ( 2 )  =  'J71' 

THE  CALL  TO  DBTAPE  REQUESTS  THE  DATA  BASE  TAPE  CONTAINING  THE 
NAMED  SATELLITE  (SAT)  AND  MODELS  (MODABBR)  AND  ASSIGNS  IT  TO 
UNIT  NUMBER  DBUNIT.  DBTAPE  SEARCHES  THE  FILE 
DBINDEX/UN=BRYANTC .  THE  WORD  NUMBERS  OF  THE  REQUESTED  MODELS 
IN  THE  DATA  RECORD  ARE  RETURNED  IN  THE  ARRAY  MRWORD  AND  THE 
FULL  MODEL  NAMES  ARE  RETURNED  IN  THE  ARRAY  MODNAME. 

CALL  DBTAPE  ( DBUNIT, SAT, NMODELS , MODABBR, MRWORD , MODNAME ) 

THE  CALL  TO  HEADIN  READS  THE  HEADER  RECORD  INTO  THE  COMMON  BLOCK 
/HEADER/.  THE  FIRST  ALTERNATE  RETURN  IS  TAKEN  IF  END-OF-FILE 
IS  ENCOUNTERED,  THE  SECOND  IF  PARITY. 

1000  CALL  HEADIN  ( *9998 , *9999  ) 

TRANS  IS  USED  TO  DO  A  BOOLEAN  TO  BOOLEAN  TRANSFER  OF  DATA  FROM 
ONE  ARRAY  TO  ANOTHER.  THIS  IS  DONE  IF  CHARACTER  DATA  IN  THE 
HEADER  IS  NEEDED. 

CALL  TRANS  ( HEADER, 40, 1 , CHEADER ) 


Figure  1.15  Sample  Program  Illustrating  PACKLIB  CALLS  (cont'd) 


1200 


1500 


9998 

i 

I 


THE  CALL  TO  DATAIN  READS  A  RECORD  OF  DATA  INTO  THE  COMMON 
BLOCK  /DATA/.  ILT,  IGMCO,  AND  IGMLT  SHOULD  BE  SET  TO  .TRUE. 
IF  LOCAL  TIME,  GEOMAGNETIC  CORDINATES,  OR  GEOMAGNETIC  LOCAL 
TIME  ARE  NEEDED,  RESPECTIVELY.  (FOR  PACKED  DATA  BASES  ONLY) 
JG  IS  THE  NUMBER  OF  DATA  POINTS  IN  THE  DATA  RECORD.  THE 
ALTERNATE  RETURNS  ARE  THE  SAME  AS  FOR  HEADIN.  IN  THIS  CASE, 
ON  END-OF-F1LE,  NORMAL  PROCESSING  RETURNS  TO  HEADIN  TO  READ 
IN  THE  NEXT  HEADER  RECORD. 

CALL  DATAIN  ( ILT, IGMCO, IGMLT, JG, *1000, *9998 ) 

CALL  TRANS  ( DATA, 40, JG, CDATA ) 

BEGINNING  OF  MAIN  PROCESSING  LOOP. 

DO  1500  J=1 , JG 

SELECT  DETERMINES  IF  DATA  POINT  J  COMES  FROM  SATELLITE  SAT. 
(APPLICABLE  TO  AE/S3-1  SATELLITES  ONLY.  FOR  SETA  SATELLITES, 
SELECT  IS  ALWAYS  TRUE.  NOTE:  SELECT  MUST  BE  DECLARED  TYPE 
LOGICAL . ) 

IF  ( . NOT . SELECT ( J, SAT ) )  GOTO  1500 


PROCESS  DATA  POINT  J. 


CONTINUE 

GO  BACK  AND  READ  THE  NEXT  DATA  RECORD. 

GOTO  1200 

END  OF  DATA  SET.  FINISH  UP  PROCESSING  (OR  GO  BACK  AND  GET 
ANOTHER  SATELLITE'S  WORTH  OF  DATA). 

CONTINUE 

STOP 

PARITY  ERROR. 


9999  STOP  'PARITY' 


Figure  1.16  PACKLIB  Routines 


PACKLIB  SUBROUTINES,  FUNCTIONS,  &  ALTERNATE  ENTRY  POINTS  WITH 
THEIR  ARGUMENT  LISTS.  (12/30/86) 


SUBROUTINES: 


SUBROUTINE 

DBTAPE 

( DBUNIT , SAT , NMODELS , MODABBR , MRWORD , MODNAME ) 

SUBROUTINE 

INIT 

( NTAPE ,  * ,  * ) 

SUBROUTINE 

HEADIN 

(*,*) 

SUBROUTINE 

HEADOUT 

(NT, HEADER) 

SUBROUTINE 

DATA IN 

( ILT, IGMCO, IGMLT, JG, * , * ) 

ENTRY 

READAT 

( ILT, IGMCO, IGMLT, JG, *, * ) 

SUBROUTINE 

DATAOUT 

(NT, DATA, JG) 

ENTRY 

WRITDAT 

(NT, DATA, JG) 

SUBROUTINE 

TRANS 

( BOOL 1 , NROWS , NCOLS , BOOL2 ) 

SUBROUTINE 

DATAMOD 

(J, CASE, SKIP) 

SUBROUTINE 

S8510RB 

SUBROUTINE 

PACKER 

( IT, IU, IA, NIA, LN, L , * , * ) 

ENTRY 

NUTP 

ENTRY 

OLTP 

(IT) 

EN'T  RY 

WRTP 

ENTRY 

PKTP 

ENTRY 

OUTP 

ENTRY 

ENTP 

ENTRY 

RETP 

ENTRY 

RDTP 

ENTRY 

UPTP 

( IT, IU, IA, NIA, LN, *, * ) 

ENTRY 

INTP 

ENTRY 

WRHEAD 

ENTRY 

RDHEAD 

( IT, IU, IA, NIA, * , * ) 

ENTRY 

BSTP 

SUBROUTINE 

REWDB 

SUBROUTINE 

DBCREAT 

( CHEADER, DBDATE ) 

SUBROUTINE 

TIMAG1 

( GL AT , GLONE , T IME , HRMNSC, GML AT , GMLONE , GMLT ) 

SUBROUTINE 

XLOCTIM 

( TIME, XLONG, HRMNSC ) 

SUBROUTINE 

HMS 

(SEC, HRMNSC) 

SUBROUTINE 

TIMAG 

( GL AT , GLONE , AM JD , GML AT , GMLONE , GMLT ) 

SUBROUTINE 

DIMAG 

( GLAT , GLONE , GML AT , GMLONG ) 

SUBROUTINE 

SOLONG 

( AM JD, ALONG) 

SUBROUTINE 

YRMODAF 

( JDATE , YEAR , MONTH , DAY ) 

SUBROUTINE 

BHWIN 

SUBROUTINE 

SKIPF 

( UNITNO, FORM, * ) 

SUBROUTINE 

COPYCF 

( FROM, TO, * ) 

SUBROUTINE 

DESPACE 

(ALPHA) 

SUBROUTINE 

DECOMMA 

(ALPHA) 

SUBROUTINE 

PACK 

(ALPHA) 
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SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

SUBROUTINE 

ENTRY 

SUBROUTINE 

SUBROUTINE 

ENTRY 

ENTRY 


LJUST  ( ALPHAI , ALPHAO , NCHAR ) 

RJUST  (ALPHAI, ALPHAO, NCHAR) 
POWER 10  (X,Y,CHSIZ, ORIENT, EXP) 
SORTS YM  (X, Y,CHSIZ, ANCHAR) 

ITOR  ( I A, A, NA ) 

RTOI  (A, I A, NA ) 

JOBINFO  (PARAM, VALUE) 

SETIND  (BOOLVAR) 

SETINF  (BOOLVAR) 

SETNINF  (BOOLVAR) 


a 


FUNCTIONS: 


LOGICAL 

FUNCTION 

SELECT 

( J,SAT) 

LOGICAL 

FUNCTION 

ALTDBF 

(SAT) 

INTEGER 

FUNCTION 

MJDATEF 

( JDATE ) 

ENTRY 

AMJDF 

( A JDATE ) 

INTEGER 

FUNCTION 

JDATEF 

( MONTH ) 

INTEGER 

FUNCTION 

MONTHF 

( JDATE ) 

INTEGER 

FUNCTION 

SEASONF 

( JDATE ) 

INTEGER 

FUNCTION 

NDTOD 

(YEAR, MONTH, DAY) 

INTEGER 

FUNCTION 

NDINYR 

( JDATE ) 

INTEGER 

FUNCTION 

NDINMO 

( JDATE ) 

REAL 

FUNCTION 

UTF 

( LT , LON ) 

REAL 

FUNCTION 

LTF 

( UT , LON ) 

REAL 

FUNCTION 

LONF 

( UT, LT ) 

INTEGER 

FUNCTION 

HHMMF 

( HOURS ) 

CHARACTER* ( * ) 

FUNCTION 

CENTER 

( ALPHA , NCHAR , F IELDWD 

CHARACTER* ( * ) 

FUNCTION 

ENCODEI 

(I,W) 

CHARACTER* ( * ) 

FUNCTION 

ENCODEF 

( X, W, D , DTRAILO ) 

CHARACTER* ( * ) 

FUNCTION 

ENCODEE 

<X,W,D) 

CHARACTER* 10 

FUNCTION 

PLOTDEV 

(  ) 

INTEGER 

FUNCTION 

LOWERF 

( LETTER ) 

REAL 

FUNCTION 

XMAXF 

(  ) 

REAL 

FUNCTION 

YMAXF 

(  ) 

LOGICAL 

FUNCTION 

IAOF 

(  ) 

LOGICAL 

FUNCTION 

IOUNITF 

(  ) 

CHARACTER* 10 

FUNCTION 

TYPE 

( BOOLVAR ) 

LOGICAL 

FUNCTION 

MINUSOF 

( BOOLVAR ) 

REAL 

FUNCTION 

SIGNUMF 

(X) 

REAL 

FUNCTION 

RGASF 

(  ) 

REAL 

FUNCTION 

DELETEF 

(  ) 

REAL 

FUNCTION 

IGNOREF 

(  ) 

1.5  Model  Development 


1.5.1  Jacchia  70  Tides 

The  early  Jacchia  models  are  based  largely  on  high  altitude 
satellite  orbital  decay  data,  and  thus  reflect  the  dominance 
of  the  diurnal  (one  oscillation  per  day)  tidal  modes  at  these 
altitudes,  ignoring  competition  from  the  semidiurnal  (twice 
per  day  oscillations)  which  occur  at  lower  altitudes.  These 
early  models  nevertheless  enjoy  the  advantage  of 
computational  simplicity  in  comparison  with  more  elaborate 
recent  models  if  only  the  total  density  is  being  computed,  as 
opposed  to  composition.  This  is  because  the  total  density 
could  be  stored  for  convenient  retrieval  from  two-dimensional 
lookup  tables.  Except  for  this  principal  weakness  of  the 
simpler  early  models  with  regard  to  their  low-altitude  local 
time  variations,  these  models  have  achieve  accuracies 
comparable  to  those  of  the  more  recent  models,  according  to 
data/model  comparisons.  The  Jacchia  70  Tides  model<712> 
represents  an  attempt  to  correct  the  local  time  variation  in 
the  Jacchia  70  model.  Polynomial  fits  in  altitude  and 
latitude  were  constructed  to  the  total  density  diurnal  and 
semidiurnal  results  of  Forbes,  et.  al.  (see  References  7  and 
12,  and  references  therein,  for  details),  and  added  to  the 
diurnally  avereged  Jacchia  70  model. 

In  statistical  evaluations  using  accelerometer  data,  as 
described  in  previous  sections  of  this  chapter,  the  Jacchia 
70  Tides  model  yielded  lower  variability  (standard 
deviations)  than  the  original  model  for  the  satellite  with 
the  most  extensive  local  time  coverage  (AE-E),  about  the  same 
variability  for  the  other  AE/S3-1  satellites,  and  somewhat 
higher  variability  for  the  SETA  series.  This  indicates  that 
probably  the  local  time  representation  in  the  Jacchia  70 
Tides  model  is  accurate,  at  least  at  low  latitudes  and  low 
solar  activity,  but  degrades  for  the  SETA  satellites  because 
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of  inaccurate  representation  of  latitudinal  and/or  solar 
cycle  dependencies  of  the  diurnal  and  semidiurnal  parameters 
A  polar-orbiting  satellite  at  high  solar  activity,  rotated  6 
hours  in  local  time  from  the  SETA  satellites,  could  possibly 
resolve  this  issue  and  lead  to  improved  modelling. 


1.5.2  Air  Force  Reference  Atmosphere  Supplements  1986 

The  Air  Force  Reference  Atmosphere  Supplements  1986  Model 
( AFRA-86 )  has  as  its  objective  the  specification  of  various 
atmospheric  properties  in  the  altitude  region  80-200  km. 
Early  on  it  was  decided  that  the  model  would  be  put  together 
from  the  MSIS-83  model  and  the  lower  altitude  model  reported 
by  Forbes  in  previous  sections  of  this  report.  Two  factors 
become  immediately  obvious: 

1)  The  Forbes  model  is  defined  only  up  to  120  km,  and 
contains  dependencies  on  only  altitude,  latitude,  and  month; 


2)  The  MSIS-83  model  is  defined  from  85  km  up  and  contains 
far  more  dependencies,  such  as  local  time,  solar  activity, 
and  geomagnetic  activity. 


Thus  it  is  apparent  that  the  Forbes  model  should  be  used  in 
the  lower  portion  of  the  AFRA-86  altitude  range  (80-200  km), 
while  the  MSIS-83  model  should  be  used  in  the  upper  portion. 


To  determine  precisely  what  the  boundaries  should  be  for 
these  respective  portions,  and  how  to  effect  the  transition 
between  them,  we  made  some  comparisons  between  the  two  models 
in  the  overlap  region  90-120  km.  In  these  comparisons  we 
deleted  from  the  MSIS-83  model  the  dependencies  not  included 
in  the  Forbes  model.  It  was  hoped  that  the  two  models  would 
be  in  close  agreement  over  a  large  enough  altitude  range  that 


the  joining  could  be  effected  by  a  small  adjustment  in  one 
or  both  models  (for  example,  modifying  selected  MSIS-83  low 
altidude  parameters)  without  compromising  either  model  in  its 
region  of  dominant  application.  Unfortunately  we  found  this 
not  to  be  the  case,  as  shown  for  example  in  Fig.  1.17  In 
general  it  was  found  that  the  Forbes-MSIS  disagreement 
worsens  as  one  goes  from  low  to  middle  latitudes,  an  effect 
apparently  related  to  differences  in  the  mid  latitude  data 
sets  used  in  their  respective  derivations.  We  therefore 
concluded  that  the  models  could  not  be  brought  into 
agreement  in  their  overlap  region  by  any  simple  adjustment 
without  jeopardizing  their  validities  in  their  separate 
altitude  regions  of  application  ( low  altitudes  for  the  Forbes 
model,  high  altitudes  for  the  MSIS-83  model). 

We  have  thus  adopted  a  model  definition  (Figure  1.18)  which 
employs  strictly  the  Forbes  model  below  a  prescribed  altitude 
hj,  the  MSIS-83  model  above  another  prescribed  altitude  h2, 
and  a  connecting  model  which  matches  continuously  and 
smoothly  the  Forbes  model  at  hj  and  the  MSIS-83  model  at  h2. 
The  value  of  h2  was  chosen  to  be  120  km,  since  it  was  found 
(Fig.  1.19)  that  the  MSIS-83  density  dependence  on  geomagneti 
activity  is  minimal  at  this  altitude,  thus  easing  the 
transition  to  the  lower  altitude  Forbes  model,  which,  as 
mentioned  previously,  contains  no  dependence  on  geomagnetic 
activity.  This  near  independence  of  the  MSIS-83  model 
density  on  geomagnetic  activity  is  evidently  due  to  the 
vanishing  of  the  coefficients  k00a  and  k20a  specifying 
the  N2  dependence  at  120  km  on  geomagnetic  activity  [see 
Hedin's  paper  [2) ,  Table  2d],  The  simplest  type  of 
connecting  function  meeting  the  requirements  of  continuous 
and  smooth  joining  at  two  points  would  be  a  cubic  polynomial 
possessing  the  required  values  and  derivatives.  Due  to 
paractical  difficulties  in  computing  derivatives,  we  have 
chosen  to  approximate  this  by  a  cubic  polynomial  matching  the 
Forbes  model  values  at  hx  -  2  km  and  hx ,  and  the  MSIS-83 
model  values  at  h2  and  h2  +  2  km. 
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Figure  1.18  Schematic  of  AFRA-86  Model 
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It  was  also  desirable,  but  not  mandatory,  that  the  connection 
model  be  consistent  with  basic  physics,  i.  e.,  satisfy 
hydrostatic  equilibrium.  The  problem  with  this  is  that  the 
solutions  to  the  equation  for  hydrostatic  equilibrium  are 
specified  by  boundary  conditions  at  only  one  altitude,  and  we 
have  to  satisfy  boundary  conditions  (the  matching  conditions 
discussed  in  the  last  paragraph)  at  two  points.  We  attempted 
to  do  this  by  modifying  the  temperature  or  molecular  weight 
so  that  the  resulting  equations  would  yield  solutions 
providing  the  required  matches. 

The  equation  for  hydrostatic  equilibrium  is: 

d( In  P )/dz  =  -Mg/(RT),  ( 1  ) 

where : 

P  =  pressure; 

z  =  altitude; 

M  =  molecular  weight; 

g  =  acceleration  of  gravity; 

R  =  universal  gas  constant 

=  8.31432xl07  g  cm2/(sec2  mole  K); 

T  =  temperature. 

The  acceleration  of  gravity  g  is  commonly  given  by: 
g  =  g8/(  l  +  z/Rp)2. 

where, 

g8  =  9.80665  m/s2; 

Rp  =  6356.766  km. 
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from  the  geometric  altitude  z  to  the  geopotential  altitude: 
y  =  (z-z1)(Rp+z1)/(Rp+z), 

where  zx  is  an  arbitrary  altitude.  Equation  (1)  then  becomes 
d( In  P)/dy  =  -Mgi/(RT), 

where  gx  is  the  acceleration  of  gravity  at  z=hx. 

It  is  convenient  to  set  zx  to  hx  and  expand  the  molecular 
weight  and  the  inverse  temperature  in  powers  of  y: 


M  =  2  nty1; 

i=0 


1/T  =  2  t^+sy^y-y^2. 


where  the  mA  and  tx  are  the  coefficients  of  the  matching 
cubic  polynomials  for  M  and  1/T,  s  is  an  adjustable 
parameter,  and  y2  is  the  geopotential  altitude  at  z  =  h2. 

Note  that  the  form  of  the  term  multiplied  by  s  preserves  the 
match  in  1/T  at  y  =  yx  =  0  and  y  =  y2.  The  solution  P  can  be 
generated  analytically  which  satisfies  the  matching 
conditions  at  both  boundaries  y*  and  y2  provided  s  is  chosen 
so  that 

s  =  [~ln(P2/P1)  -  (g1F/R)]/(g1G/R), 


where : 


Px  =  pressure  at  yx; 


P2  =  pressure  at  y2; 


6 

f  =  y2  *  Q1Y2  ; 

1=0 

3 

g  =  y25  ^  "'xriy2i, 

1=0 

where: 

qi  =  (  2  nijti.j )/(  i  +  1 ) ; 

j=0 

Ti  =  l/(  1  +  5 )-2/( i+4 )  +  l/( 1+3 ) . 

Given  P,  T,  and  M  over  the  interval  to  y2  we  can 
calculate  the  density  p  from  the  perfect  gas  law: 

P  =  PM/(RT). 

A  similar  formulation  can  be  developed  for  adjusting  the 
molecular  weight  M  rather  than  1/T,  i.  e.,  we  add  a  term  to 
the  polynomial  for  M  containing  a  coefficient  which  is 
adjusted  to  make  the  required  match. 

We  examined  these  approaches  for  various  cases.  Unfortu¬ 
nately  we  found  some  cases,  as  in  Figure  1.20,  where  either 
the  temperature  had  to  go  through  at  minimum  at  an  unrealis¬ 
tically  high  altitude  (normally  it  should  minimize  near  90 
km)  or  the  mass  attained  unrealistically  high  or  low  values. 
Broadening  the  transition  region  did  not  help,  nor  was  it 
felt  feasible  to  add  more  terms  to  the  expansions,  since  this 
would  complicate  the  calculations  unduly  and  most  likely 
would  simply  result  in  solutions  containing  more  severe 
undulations. 


46 


Figure  1.20  Hydrostatic  Equilibrium  Model  in  Transition 
Region  90-120  km 


Therefore  we  have  chosen  to  drop  the  requirement  of  static 
equilibrium  and  define  the  Forbes-MSIS  transition  strictly  by 
matching  cubic  polynomials.  After  examination  of  typical 
results  we  chose  104  km  as  the  lower  boundary  of  the 
transition  region.  This  amounts  to  a  compromise  between  110 
km  suggested  by  Forbes  as  the  upper  limit  of  validity  of  his 
model  and  a  more  conservative  upper  limit  at  100  km,  which 
arose  out  of  our  preliminary  studies  (we  found  in  particular 
severe  disagreement  between  Forbes  and  the  U.  S.  Standard  76 
Atmosphere  at  110  km,  45  degrees  latitude,  as  shown  in 
Figure  1.21). 

In  conclusion,  we  we  give  the  final  definition  for  the  Air 
Force  Reference  Atmosphere  1986  (AFRA-86): 

Region  1  (70-104  km):  Forbes  model; 

Region  2  (104-120  km): 

compute  the  density,  temperature  and  molecular 
weight  from  cubic  polynomials  matching  the  Forbes 
model  at  102  km  and  104  km,  and  matching  the 
MSIS-83  model  at  120  km  and  122  km;  then  compute 
the  pressure  from  the  perfect  gas  law; 

Region  3  (120-200  km):  MSIS-83  model. 


Study  of  Statistical  Techniques 


Applied  to  Atmospheric  Density  Data 

1.6.1  Introduction 

It  was  suggested  that  the  neutral  mass  density  to  model 
ratios  may  be  corrrelated  in  time.  An  investigation  was 
undertaken  to  evaluate  this  possibility  and  to  assess  it's 
effect  on  the  statistics  currently  being  calculated  for  Air 
Force  Reference  Atmosphere  purposes.  Two  approaches  were 
used:  one  was  to  do  a  one  way  random  effects  analysis  of 
variance,  separating  the  variation  into  between  and  within 
orbit  components,  the  second  was  to  calculate  the  auto¬ 
correlation  coefficients  for  various  lags  and  compare  them 
with  their  standard  errors.  The  Durbin-Watson  statistic  for 
serial  correlation  was  also  calculated. 


1.6.2  One  Way  Random  Effects  Analysis  of  Variance 

The  STAT  program  was  modified  to  do  a  one  way  random  effects 
analysis  of  variance  calculation.  The  Jacchia  1971  model 
ratios  were  analyzed  from  the  SETA-1  data  set  spanning  days 
79  through  100  of  1979.  The  data  were  binned  on  4  Kp  bins  (0 
to  1+,  2-  to  3+,  4-  to  5+,  and  6-  to  9)  and  5  geographic 
latitude  bins.  The  variation  within  each  bin  was  divided 
into  between  and  within  orbit  variance  components. 


The  one  way  random  effects  model  is: 


Yij  =  M  +  «i  +  €ij 

where  Yij  is  the  mean  ratio, 

u  is  the  overall  mean  in  each  bin 

oq  is  a  random  component  due  to  orbit  to  orbit 
variability  and, 

cjj  is  a  random  component  due  to  variability 
within  orbits 

(See  References  15  and  16) 

Note  that 

E[<Xi]  =  0,  V[ocJ  =  a*2 

E^]  =  0,  VCe.j]  =  a2 

so  that 

EtYij]  =  M 
VtYij]  =  ^(x2  ♦  °-2 

where  E [ . ]  is  the  expected  value,  and  V[.]  is  the 

variance.  We  want  to  test  if  cft2  =  0.  If  the  data  are 
correlated  in  time,  one  would  expect  c2  to  be  much  less  than 


We  construct  the  following  analysis  of  variance  table. 


One  way  Random 

Effects 

Analysis  of 

Variance  Table 

Source 

Sum  of 
Squares 

Degrees 

Freedom 

Mean 

Square 

Expected 

Mean  Square 

Total 

SST 

N-l 

Between 

SSB 

MSB 

a2  + 

nea0<2 

Within 

SSW 

N-n0 

MSW 

a2 

where 

N  =  total  number 

of  points. 

nQ  =  number  of  orbits. 

MSB  = 

SSB 

and. 

MSW  = 

SSW 

• 

N-n0 

The  total 

sum  of  squares  can 

be  partioned  into 

its  between 

and  within 

components . 

SST  = 

ni 

*1.1  ( Yij 

-  y)2 

-  2n° 

nj  <^J  * 

r>2  ♦  Zj. 

n°  ^  nj/ 
i  2i-i  (Yu  - 

y.S 

=  SSB 

+  SSW 

where  y.j  is  the  average  over  i  =  l,  nj . 

Note  that  E[MSB]  =  a2  +  n,,^2 
and  E[MSW]  =  <r2. 

n„  is  equal  to  or  less  than  the  number  of  points  per  bin  per 
orbit  with  equality  holding  when  there  are  an  equal  number  of 
points  per  bin  per  orbit  (See  Reference  17). 

To  test  that  cr^2  =  0,  we  can  test  the  F  ratio, 

F  =  MSB 
MSW 

for  significance.  Table  1.6  shows  the  calculation  for  the 
SETA-1,  Jacchia  71  data  set.  The  value  in  the  column  "l-P" 
is  the  probability  of  exceeding  the  calculated  F  ratio. 


Table  1.6:  SETA-1  Jacchia  71 
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Table  1 . 6  cont 


It  is  readily  apparent  that  nearly  every  case  showed 
statistical  significance  leading  one  to  reject  the  hypothesis 
that  aw2  =  0.  The  tabulated  values  of  S.D.B  and  S.D.W  are 
the  square  roots  of  the  variance  components  cra2  and  a2  (i.e. 
standard  deviations)  respectively.  It  can  be  seen  that  the 
standard  deviation  within  orbits  is  less  than  the  standard 
deviation  between  orbits  which  would  lead  one  to  conclude 
that  the  data  are  correlated  within  orbits. 


1.6.3  Calculation  of  the  Autocorrelation  Coefficients 

Next,  20  random  samples  of  100  consecutive  data  points  were 
taken  from  the  SETA-1,  Jacchia  71,  21  day  data  set.  The 
Durbin-Watson  statistic  was  calculated  and  the  auto¬ 
correlations  and  their  standard  errors  were  calculated  for 
lags  up  to  lag  25  for  each  sample. 

The  Durbin-Watson  statistic  is  defined  as: 

*i-2  (Yi  "  Yi-i) 

a  =  - 

(Yi  -  r>2 

where  yA  is  the  model  ratio. 

In  this  case  we  are  testing  for  p  >  0  and  small  values  of  d 
are  significant  (See  Reference  14). 

The  sample  auto  correlation  coefficient  at  lag  1  is  defined 
as 

^i1.!  (Yi  ~  ~Y )  (Yi.i  "  Y) 

ri  - 

^i-i  (Yi  -  y}2 


where  yt  is  the  model  ratio. 


The  large  lag  standard  error  of  the  auto  correlation 
coefficient  is  defined  as: 


V[rJ  =  -  (1  +  21?ml  r,2  )  l>q 

N 

where  r*  is  the  sample  auto  correlation  coefficient  at  lag  i. 
(see  Reference  13) 

The  results  of  this  calculation  are  given  in  Table  1.7.  The 
Durbin-Watson  statistic  was  significant  at  the  5%  and  1% 
levels  for  all  20  samples  indicating  evidence  of  positive 
serial  correlation.  Inspecting  Table  1.7,  one  notices  that 
the  auto  correlations  die  out  so  that  they  are  typically  on 
the  order  of  one  or  two  times  their  standard  errors  by  fourth 
or  fifth  lag.  This  would  indicate  that  to  achieve  indepen¬ 
dence  in  sampling  one  should  sample  at  most  only  every  fourth 
or  fifth  data  point. 


1.6.4  Effect  of  Autocorrelation  on  Analysis  of  Satellite  Data 


Assume  that  on  average  N  samples  are  taken  during  each  orbit, 
and  that  nQ  orbits  of  data  are  collected.  Assume  the  n 
samples  per  orbit  are  "perfectly”  correlated.  Since  the  n 
samples  per  orbit  are  essentially  repeats  we  obtain  a 
standard  deviation  of 


nQ  _  _ 2 

n  Sj-i  (y.j  -  y) 

N  -  1 


where  N  =  n.n„ 
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Table  1.7 


Autocorrelation 
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-«vV: 
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if  we  restrict  ourselves  to  one  point  sampled  at  random  per 
orbit  we  will  find  instead  a  standard  deviation  of 


s/  = 


Sjti  ( Yj  -  y): 

n„-l 


The  relative  standard  deviations  in  the  two  cases 
1)  oversampling,  and  2)  actual,  are  thus  in  the  ratio 


Take  a  typical  case  of  30  orbits  with  4  samples  per  orbit 
The  apparent  standard  deviation  to  the  true  standard 
deviations  are  in  the  ratio 


or  1.3%  lower  than  actual. 


1.6.5  Conclusion 

We  conclude  that  while  the  data  do  appear  to  be  correlated  in 
time,  this  has  little  or  no  impact  on  the  magnitude  of  the 
standard  deviations  as  they  are  currently  being  calculated. 
However,  it  does  inpact  the  "degrees  of  freedom"  associated 
with  the  standard  deviations  which  will  impact  any  confidence 
or  prediction  intervals  that  may  be  calculated  from  these 
statistics. 
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1.7  AE-E  Local  Time  Effect 


Since  the  AE/S3-1  data  is  characterized  by  greater  local  time 
coverage  than  the  SETA-1  data,  it  has  been  suspected  that 
unmodeled  local  time  variations  in  the  density  contribute 
significantly  to  the  higher  variabilities  (standard 
deviations)  of  the  AE/S3-1  measured  to  model  ratios. 

Therefore  an  analysis  of  variance,  similar  to  that  described 
in  Section  1.6,  was  conducted  for  the  AE-E  data,  which  has 
the  most  extensive  local  time  coverage,  to  assess  the 
significance  of  unmodeled  local  time  variations,  and  to  see 
if  elimination  of  these  variations  could  bring  the  AE/S3-1 
and  SETA  variabilities  into  agreement.  Results  are  given  in 
Tables  1.8,  1.9  for  the  data  to  Jacchia  71  model  ratio  and 
the  data  to  MSIS  83B  model  ratio.  These  results  show  that, 
while  the  unmodeled  local  time  variations  are  significant  for 
the  AE-E  data,  their  elimination  reduces  the  AE-E  variability 
only  to  around  12%,  still  well  above  the  SETA  variabilities. 
Note  also  that  the  unmodeled  local  time  variations  are  more 
significant  for  the  Jacchia  71  model  ratios  than  for  the  MSIS 
83B  ratios,  indicating  that  the  local  time  variations  for  the 
AE-E  data  are  better  represented  by  the  MSIS  83B  model  than 
by  the  Jacchia  71  model. 
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Table:  1.8  ONE 

WAY  ANALYSIS 

OP  VARIANCE 

FOR  AE-E 

LOCAL  TIME  EFFECT 

JACCHIA 

71  MODEL 

MODEL: 

LN(  R(  I ,  J ) )  =  MU  +  LT  ( I ) 

+  E(I,J) 

1-1,24,  J-1,N(I) 

SOURCE  OP 

SOM  OF 

DEGREES  OF 

MEAN 

F  PR(F>1009 . 8) 

VARIATION 

SQUARES 

FREEDOM 

SQUARE 

STATISTIC 

TOTAL 

826.262 

35628 

.02306 

LOCAL  TIME 

326.112 

23 

14.178 

1009.8  0.0* 

RESIDUAL 

500.150 

35605 

.01404 

*  WITHIN  OUR  ABILITY  TO  COMPUTE  THE  RELEVANT  PROBABILITY  LEVEL. 

OVERALL  MEAN  MU  =  .0144 

TOTAL  VARIANCE  =  .02306 

TOTAL  STANDARD  DEVIATION  -  .1523 
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I 

LOCAL  TIME 

N(I) 

MU  + 

VARIANCE 

S.D. 

(HOURS) 

LTd) 

WITHIN 

WITHIN 

1 

0-1 

1138 

.1501 

.01861 

.1364 

2 

1-2 

1111 

.1378 

.02030 

.1425 

3 

2-3 

1116 

.1417 

.01552 

.1246 

4 

3-4 

1138 

.1772 

.01859 

.1363 

5 

4-5 

1061 

.1754 

.01787 

.1336 

6 

5-6 

1662 

.1919 

.01773 

.1332 

7 

6-7 

1980 

.1689 

.01357 

.1165 

8 

7-8 

2183 

.1647 

.01356 

.1164 

9 

8-9 

1996 

.1149 

.01142 

.1069 

10 

9-10 

2018 

.0433 

.01182 

.1087 

11 

10-11 

1797 

.0090 

.01130 

.1063 

12 

11-12 

1587 

-.0239 

.00897 

.0947 

13 

12-13 

1593 

-.0432 

.00856 

.0926 

14 

13-14 

1485 

-.0764 

.01107 

.1052 

15 

14-15 

1449 

-.1029 

.01215 

.1102 

16 

15-16 

1372 

-.0698 

.00983 

.0992 

17 

16-17 

1320 

-.0559 

.01327 

.1152 

18 

17-18 

1338 

-.0602 

.01205 

.1098 

19 

18-19 

1401 

-.0143 

.01845 

.1358 

20 

19-20 

1452 

.0778 

.01506 

.1227 

21 

20-21 

1389 

.1322 

.01183 

.1088 

22 

21-22 

1397 

.1252 

.02194 

.1481 

23 

22-23 

1397 

.1251 

.01643 

.1282 

24 

23-24 

1249 

.1224 

.01564 

.1251 

RESIDUAL 

VARIANCE 

-  .01404 

RESIDUAL 

STANDARD 

DEVIATION 

»  .1185 

,%  •> 
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Table:  1 .9  ONE  WAY  ANALYSIS  OP  VARIANCE  FOR  AE-E  LOCAL  TIME  EFFECT 


MSIS  83B  MODEL 


MODEL:  LN{R(I,J))  =  MD  +  LT(I)  +  E(I,J)  1-1,24,  J-1,N(I) 


SOURCE  OF  SUM  OF  DEGREES  OF  MEAN  F  PR (F> 3 19. 6) 

VARIATION  SQUARES  FREEDOM  SQUARE  STATISTIC 


TOTAL  619.688  35628  .01739 

LOCAL  TIME  106.082  23  4.612  319.6  0.0* 

RESIDUAL  513.606  35605  .01443 


WITHIN  OUR  ABILITY  TO  COMPOTE  THE  RELEVANT  PROBABILITY  LEVEL. 


OVERALL  MEAN  MU  -  .0144 

TOTAL  VARIANCE  -  .01739 

TOTAL  STANDARD  DEVIATION  -  .1319 


I  LOCAL  TIME 

N(I) 

MO  + 

VARIANCE 

S.D. 

(HOURS) 

LT(I) 

WITHIN 

WITHIN 

1 

0-1 

1138 

.0554 

.02179 

.1476 

2 

1-2 

1111 

.0447 

.02195 

.1481 

3 

2-3 

1116 

.0480 

.01319 

.1148 

4 

3-4 

1138 

.0747 

.01650 

.1284 

5 

4-5 

1061 

.0433 

.02015 

.1419 

6 

5-6 

1662 

-.0027 

.01666 

.1291 

7 

6-7 

1980 

-.0347 

.01594 

.1262 

8 

7-8 

2183 

-.0362 

.01309 

.1144 

9 

8-9 

1996 

-.0461 

.01143 

.1068 

10 

9-10 

2018 

-.0584 

.01391 

.1144 

11 

10-11 

1797 

-.0316 

.01320 

.1149 

12 

11-12 

1587 

-.0258 

.00992 

.0996 

13 

12-13 

1593 

-.0348 

.00798 

.0893 

14 

13-14 

1485 

-.0363 

.01213 

.1101 

15 

14-15 

1449 

-.0240 

.01165 

.1079 

16 

15-16 

1372 

.0416 

.00798 

.0893 

17 

16-17 

1320 

.0467 

.01222 

.1105 

18 

17-18 

1338 

.0312 

.01118 

.1057 

19 

18-19 

1401 

.0569 

.01589 

.1260 

20 

19-20 

1452 

.1235 

.01442 

.1200 

21 

20-21 

1389 

.1399 

.01256 

.1121 

22 

21-22 

1397 

.0795 

.02224 

.1491 

23 

22-23 

1397 

.0421 

.01868 

.1366 

24 

23-24 

1249 

.0288 

.01933 

.1390 

RESIDUAL 

VARIANCE 

-  .01443 

RESIDUAL 

STANDARD 

DEVIATION 

-  .1201 
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2 . 0  Ionospheric  Scintillation  Data  Processing  Systems 

2 . 1  Overview 

Two  radio  wave  scintillation  data  processing  systems, 
developed  by  SRI  International,  known  as  AFSATCOM  (SDS)(1'2) 
and  Fleetsatcom  ( FLT ) 131  routinely  process  radio  wave 
scintillation  data  from  near  stationary  polar  beacons  and 
geostationary  satellites,  respectively.  Data  are  collected 
using  an  Intel  8080  microprocessor  system  at  Goose  Bay, 
Labrador,  Thule,  Greenland,  Tromso,  Norway,  and  Malvik, 
Scotland,  for  SDS  and  at  Ascension  Island  for  the  FLT  system. 
The  format  of  the  field  tapes  is  shown  in  Table  2.1.  Data 
reduction  and  analysis  takes  place  at  AFGL  and  involves 
unpacking  the  raw  data,  separating  it  into  its  phase  and 
intensity  components,  extracting  that  portion  of  total  phase 
ascribable  to  propagation  effects,  and  deriving  certain 
statistical  properties  of  the  propagation  data  (See  Table  2.2 
and  Figures  2.1,  2.2,  2.4,  and  2.5).  In  the  process,  the 
data  are  spectrum  analyzed  and  plots  of  the  power  spectrum  of 
the  data  are  part  of  the  output  of  the  SDS  and  FLT  systems. 
Statistics  on  each  block  of  data  are  accumulated,  usually  on 
a  monthly  basis,  and  summary  statistics  are  generated  using 
the  STATSDS  and  STATFI.T  programs  (See  Figures  2.3  and  2.6). 

In  addition,  SUATEK  plots  of  parameters  of  interest  are  part 
of  the  period  end  processing  (See  Table  2.3). 


i  m  y  *, 


i 


Byte  # 


Meaninc 


11-12 


16-27 


31-46 


48-49 

50-51 


53-63 

64-1663 


record  type  (1  =  data,  2  =  calibration) 
pass  number 

record  serial  number  (l.s.  byte  first) 

seconds 

minutes 

hours 

days 

year 

system  configuration  word 
receiver  status  (l.s.  byte  first) 
mode 

chart  speed 
chart  flag 

most  recent  HPIB  tuning  command  string 
(ASCII) 

number  of  channels  in  A/D  list, 
number  of  samples/burst 
offset  into  list  for  cosine  channel 
offset  into  list  for  sine  channel 
list  of  A/D  channels  in  sequence 
divide  factor  for  10  ms  clock 
hour,  minute  to  open  window 

hour,  minute  to  close  window 

hour,  minute  to  close  window 

blank,  spare  space 

A/D  converter  samples,  2  bytes/sample 


(INTEL  8080  format:  16  bit,  two's  complement, 

least  significant  byte  first) 


Table  2.2  Unpacked  Data  Format 


Word 

Meaning 

1-800 

(intensity,  phase)  pairs* 

801 

time  in  seconds 

802 

day  number 

803 

sampling  frequency,  Hz 

804 

channel  #1  frequency,  Hz 

805 

channel  #2  frequency,  Hz 

806 

channel  #3  frequency,  Hz 

807 

pass  number 

808 

record  serial  number 

809 

system  configuration  word 

810 

receiver  status 

811 

mode 

*  Note: 


For  mode  3  data,  the  intensity  (I)  and  phase  pairs 
for  the  three  channels  are  interleaved 
(i.e.  (In/  Pli)<  (l3i»^3i)>  i  =  l,133), 

with  two  extra  "dummy"  values  at  the  end  to  pad 
the  array  out  to  800  words. 


Note:  Mode  5  data  consists  of  intensity  data  only. 
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Table  2.3  SDS  and  FLT  Plots 


i 


r  - 


K 


K 


y. 


y. 


K 


i£ 


K 


V 

■A 


SDS  Plots 


1. 

Kp 

vs 

Magnetic 

Local 

Time 

2. 

RMS  Phase  Deviation 

vs 

Magnetic 

Local 

Time 

3. 

Phase 

Spectral 

Slope 

vs 

Magnetic 

Local 

Time 

4. 

Phase 

Spectral 

Strength 

vs 

Magnetic 

Local 

Time 

5. 

s4 

vs 

Magnetic 

Local 

Time 

6. 

Intensity  Spectral  Strength 

vs 

Magnetic 

Local 

Time 

7. 

Log^  ( 

Decorrelation  Time) 

vs 

Magnetic 

Local 

Time 

( Night 

and  Day ) 

8. 

Phase 

Spectral 

Slope 

vs 

s4 

9. 

Phase 

Spectral 

Slope 

vs 

RMS  Phase  Deviation 

10. 

Phase 

Spectral 

Slope 

vs 

Phase  Spectral 

Strength 

FLT  Plots 


1. 

Decorrelation 

Time 

vs 

S4  ( Pre  and  Post  Midnight ) 

2. 

Phase  Spectral 

Slope 

vs 

RMS  Phase  Deviation 

3. 

Decorrelation 

Time 

vs 

Phase  Spectral  Strength 

4. 

Phase  Spectral 

Slope 

vs 

Phase  Spectral  Strength 

83 


Figure  2.2  SDS  Data  Reduction  System 


Figure  2.3  SDS  Period-end  Statistics  Processing 


N 


Figure  2.4  FLT  Ionospheric  Scintillation  Data  Processing  System 


Figure  2.6  FLT  Period-End  Statistics  Processing 


2.2  Description  of  Procedures 


Procedures  were  written  to  perform  all  of  the  SDS  and  FLT 
data  processing  functions  (plus  some  additional  functions  not 
previously  available,  useful  in  diagnosing  anomalous  data). 
Each  of  the  processing  functions  is  now  operationally  much 
simpler,  and  any  maintenance  and/or  program  changes  are 
transparent  to  the  user. 


2.3  Discrete  Fourier  Transform  Window/MEM  Study 

The  process  of  Fourier  analysis  is  subject  to  a  number  of 
contaminating  effects.  To  insure  reliability  of  conclusions 
from  the  resulting  spectra,  it  is  important  to  minimize  the 
impact  of  such  spurious  effects. 

Leakage  is  one  of  the  major  contributors  to  contamination  of 
Fourier  spectra.  It  is  a  consequence  of  the  fact  that,  in 
practice,  one  deals  with  data  samples  of  finite  length  which 
the  Fourier  processor  treats  as  being  periodic.  The  computed 
spectrum  then  corresponds  to  that  of  a  finite  duration  sample 
repeated  periodically  from  -«•  to  +».  The  cause  of  the 
leakage  effect  is  the  convolution  of  the  rectangular  "boxcar" 
function  representing  the  finite  time  duration  sample  with  an 
assumed  signal  extending  from  -w  to  +•».  In  the  frequency 
domain,  the  consequence  of  this  effect  is  that  spectral  power 
at  one  frequency  contributes  a  certain  amount  of  spurious 
power  at  other  frequencies  in  addition  to  whatever  power 
already  exists  at  those  frequencies.  Depending  on  the 
relative  magnitudes  of  the  already  occurring  and  spurious 
power  at  a  given  frequency,  the  computed  spectrum  may  be 
seriously  in  error. 


Leakage  has  its  most  serious  impact  in  those  cases  in  which 
it  is  necessary  to  observe  weak  contributions  at  a  certain 
set  of  frequencies  in  the  presence  of  relatively  high  level 
contributions  at  nearby  frequencies.  The  contribution  of 
spurious  power  at  one  frequency  to  its  neighbor  goes  as 
sin(x)/x,  where  x  is  proportional  to  the  difference  in 
frequencies.  The  1/x  fall-off  implies  that  nearby  neighbors 
are  most  severly  affected.  In  the  current  application,  the 
most  serious  contamination  occurs  in  the  case  of  very  steep 
roll-offs  at  relatively  low  frequencies. 

The  standard  approach  to  the  problem  of  leakage  is  the  use  of 
window  functions. <4)  Their  utility  lies  in  their 
suppressing  the  magnitude  of  the  artifact  discussed  above. 
However,  a  window  is  itself  an  artifact  and  has  its  own 
disadvantage,  i.e.  the  transform  of  the  window  function  is 
convolved  with  the  unwindowed  spectrum,  tending  to  smear 
and/or  broaden  individual  features  of  the  resulting  spectrum. 
The  result  is  a  loss  of  spectral  resolution.  In  general,  the 
better  a  window  function  suppreses  leakage,  the  poorer  the 
resulting  resolution.  Hence,  a  tradeoff  is  necessary 
depending  on  the  character  of  the  data  being  analyzed. 

Another,  approach  is  to  replace  discrete  Fourier  analysis, 
with  its  Implicit  assumption  of  periodicity,  by  the  Maximum 
Entropy  Method  (MEM).<5-6)  The  MEM  avoids  the  assumption  of 
periodicity  of  the  Fourier  method. 

To  assess  the  effect  of  different  windows  on  the  current 
application  and  to  make  a  comparison  of  the  MEM  with  the 
windowed  DFT,  a  synthetic  signal  and  Tromso  Day  61  were 
processed  by  a  suite  of  windows  and  the  MEM.  The  candidate 
windows  are  given  in  Table  2.4.  As  a  result  of  this  study, 
and  at  the  request  of  the  initator,  the  -74  dB  4  term 
Blackman-Harris  window  was  implemented  in  the  SDS  system  as 
an  option. 


Table  2.4  DFT  Window  Candidates 


1 .  No  window 

2.  End  matching 

3.  cos 

4.  cos2 

5 .  cos3 

6.  cos4 

7.  -92  dB  4  term  Blackman-Harris 

8.  -74  dB  4  term  Blackman-Harris 

9.  4  sample  Kaiser-Bessel 

10.  Subtract  mean 


2.4  Statistics  Programs 


There  are  two  statistics  programs  wich  produce  period  end 
reports  from  the  merged  statistics  files,  namely  STATSDS  and 
STATFLT  (See  Figures  2.3  and  2.6).  STATFLT  is  an  entirely 
new  program,  while  STATSDS  was  rewritten  to  simplify  program 
maintenance.  Both  programs  produce  a  listing  of  the 
accumulated  statistics  on  the  merged  statistics  files  and  a 
report  summarizing  the  accululated  statistics  for  the  period 
of  interest.  A  sample  output  of  the  STATSDS  program  is  given 
in  Table  2.5.  STATFLT  performs  a  similar  function  to 
STATSDS,  but  for  the  FLT  system.  Since  there  are  two 
separate  merged  statistics  files  in  this  case,  one  for 
intensity  and  one  for  phase,  STATFLT  sorts  the  files  by  day 
number  and  start  time  and  checks  for  duplicate  and/or  missing 
intensity  and/or  phase  records.  The  output  of  STATFLT  is 
similar  to  that  of  STATSDS. 


Table  2.5  Sample  STATSDS  Output 


Table  2.5  (Cont.)  Sample  STATSDS  Output 
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2.5  SDS  Scintillation  Processing  System  User's  Guide 

The  SDS  Scintillation  Processing  System,  originally  developed 
by  Stanford  Research  Institute  International,  is  used  in 
order  to  process  radio  wave  scintillation  data  collected 
using  transmissions  from  beacons  on  board  satellites. 
Procedures  that  are  needed  to  acquire  the  field  tapes,  and 
the  subsequent  processing  for  plots,  statistics,  and 
evaluation  are  described. 

2.5.1  Introduction 

It  is  the  intent  of  this  document  to  describe  the  input 
necessary  to  run  the  various  processing  capabilities  of  the 
SDS  Scintillation  Processing  System.  It  is  not  the  intent 
of  this  document  to  describe  in  detail  the  mathematical 
computations  involved.  For  those  details,  the  reader  is 
referred  to  Reference  2. 

All  of  the  SDS  processing  functions  are  now  handled  by  NOS 
procedures.  Hopefully,  this  will  provide  for  ease  of  use  by 
the  user.  A  description  of  the  input  necessary  to  run  the 
SDS  Scintillation  Processing  System  follows. 


2.5.2  SDS  Tape  Scan 

Field  tapes  are  submitted  first  to  the  TPSCAN  procedure, 
which  provides  a  printed  listing  of  the  contents  of  the  files 
on  that  tape.  Figure  2.7  is  a  typical  tape  scan  output.  The 
start  and  end  times  for  one  file  have  been  identified.  The 
TPSCAN  program  also  writes  a  file,  TAPE19,  which  may  be  SAVED 
for  use  in  later  processing.  The  TAPE19  file  contains 
information  used  in  removing  an  unwanted  receiver  artifact 


wyim  ir*  ir«ir»  c-«u»  ** 


1 


from  the  intensity  and  phase  spectra  (See  Reference  7). 

Input  to  TPSCAN  consists  of  the  visual  serial  number  (VSN)  of 
the  tape . 

Input  to  TPSCAN: 

Card  Columns 

1  1-6 

A  sample  deck  setup  is  shown  in  Figure  2.8.  This  deck  is 
available  as  CCSDS1/UN=BRYANTW. 


Description 
VSN  of  field  tape 


•<T0P  OF  FILE) 

BRYANT, WOO.  SDS  TAfC  SCAN 
USCR, BRYANT V,BRYANTU. 

CHAR6C,550t ,4643. 

6£T,fR0Crit/UN«BRXANT«. 

LABCL , TAPE t ,VSN*M9l ,NT.O=Fe.F=S,LB*KU,PO=R. 
BEGIN.  TPSCAN, ,T 194044. 

EXIT. 

EOR 

TR4044 

<B0TT0N  OF  FILO 


Figure  2.8:  An  Exaaple*  of  Deck  Setup  for  Tape  Scan 
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2.5.3  Normal  SDS  Processing 


A  complete  description  of  the  data  processing  involved  in 
normal  SDS  processing  is  given  in  Reference  2.  Briefly,  the 
data  are  read  from  a  field  tape,  are  unpacked  and  reduced. 

The  reduced  data  is  spectrum  analyzed  and  various  statistics 
are  calculated.  The  output  of  the  normal  SDS  processing 
procedure  consists  of  a  printed  listing  of  the  statistics  for 
the  files  and  time  periods  of  interest.  Plots  of  intensity 
and  phase  spectra  are  produced.  Either  micro-fiche  or 
penplots  are  available.  Finally,  the  calculated  statistics 
are  written  to  TAPE3  which  may  be  SAVEd  for  later  processing 
(See  Sections  2.5.4  and  2.5.5). 


Input  for  Normal  SDS  Processing: 


Card  Description 

1  VSN  of  field  tape 

2  file,  station  #,  channel,  year,  window  option. 

3  start  time  (day  number,  hour,  minute,  second) 

4  end  time  (day  number,  hour,  minute,  second) 

The  file  is  the  actual  file  number  from  a  tape  scan  run. 

The  station  numbers  are  as  follows: 

Station  #  Location 


2 

3 

4 

5 


Goose  Bay 
Thule 
Tromso 
Malvik 


Set  the  window  option=0  for  no  windowing,  or  window  option=l 
for  the  -74  dB  4-term  Blackman-Harris  window  (See  Reference 
4).  All  quantities  are  INTEGER  and  in  free  format  (i.e. 
list  directed  format),  except  for  VSN  which  is  alphanumeric. 


Input  for  multiple  file  processing  should  be  separated  by  end 
of  record  marks  (EORs).  The  default  medium  for  SDS  plots  is 
microfiche.  Pen  plots  may  be  obtained  by  specifying  PEN=YES 
on  the  BEGIN, SDS  statement  (i.e.  BEGIN, SDS, , PEN=YES ) .  A 
sample  deck  setup  is  given  in  Figure  2.9.  This  deck  is 
available  as  CCSDS2/UN=BRYANTW. 


2.5.4  Merging  SDS  Data  Files 

The  TAPE3  data  files  saved  during  normal  SDS  processing  may 
be  merged  using  the  MERGIT  procedure.  The  easiest  way  to  do 
this  is  to  run  the  MERGIT  procedure  interactively.  An 
example  of  this  is  given  in  Figure  2.10.  The  user's  typeins 
have  been  underlined  in  Figure  2.10.  Note  that  the  user  must 
GET  the  procedure  file  PROCFIL/UN=BRYANTW  before  executing 
the  procedure  MERGIT. 


2.5.5  End  of  the  Month  Statistics 

End  of  the  month  statistics  on  the  merged  data  files  are 
produced  by  the  STATS  procedure.  A  listing  of  the  data  on 
the  merged  file  is  produced  along  with  the  report  shown  in 
Figure  2.11.  A  plot  file  suitable  for  plotting  by  SUATEK  is 
also  produced.  This  plot  file  contains  geomagnetic  local 
time,  Kp,  rms  phase  deviation,  phase  spectral  slope,  phase 
spectral  strength,  S4,  intensity  spectral  slope,  intensity 
spectral  strength,  decor-  relation  time,  logl0( decorrelation 
time),  and  a  day/night  indicator.  Day  is  considered  to  be 
geomagnetic  local  time  between  the  hours  of  0900  to  2100. 


HE  KM  KtfSgfglOTOSPlPSWWVWl1  W  ^  WtHWf  BBWWWWIffWMfJ*  \y*^  'JhJ'fote  rjnanvjw  >jmirmr*.>r*  tr>«  -.«  - 


•<TOf  Of*  fIL€> 

e«{VAf<T .1*1 500 .  NORftAL  SOS  PROCESSING 
»S£R.RRYANTV,BKYANTV. 

CHARGE .5501 ,4643. 

SET  ,PRQCFIl/UNsBRYANTU. 

SET .TAPE1 9*T1 94Q44/UN*&RYAHTW. 

LABEL,  TAPE  7,  V$N*fl91 ,NT,DsRE,f*S,CB»KU,RO*R. 
BE6IN.SBS. 

REPLACE. TAP£3»T4044F1 . 

BEGIN. SDS. 

REPLACE. TAPE3*T4044f4. 

BEGIN.SOS . 

REPLACE.  TAPC3*T4044F7. 

EXIT. 

COR 

TR404  4 
1.4.1  .f;4 
61 .21,3'  .0 
42.1.  tO.  44 
COR 

1R49-M 
*,4,1.04 
42.9.15.52 
62, 1 .5. 1 4,46 
COR 

TR4044 
7.4,1 .04 
62,21.21,13 
63,1,10,49 
< BOTTOM  Of  FILE) 
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Figure  2.9:  An  Example  of  Deck  Setup  fop  Normal  SDS  Processing 
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<wt.orocfil/uwbrvBwtti 

WELCOME  TO  THE  SIS  MERC1N6  PROCEDURE. 

CnUr  PFN2  SOS  PILE  NAME  TO  *€  HERBERT  1404491 
Eittr  PFN3  SIS  OR  OCR  MERGES  FILE  NAME  T  14Q44T4 
E»ltr  PFN4  NEW  NER6EI  FILE  MAMET  14044 

15. 17. 31.  USER  BATFILE  PROCESSED. 

15.17.31.  COMMENT .  mmmwmwwmmmmm 

15. 17. 31.  COMMENT.  4  SOS  NERSINS  PROCEDURE  • 
15.17.31  .COMMENT. 

15.17.31 . BE6IM»SETPfN, .T4044F1 . 
13.t?.33.IF,iO«.EQ.m,LAB£U. 

15. 17. 33.  BET .T4044F1/NA. 

1 5. 17. 33.  ELSE, LABEL 1 . 

15. 17.33.  ENDIF , LABEL  1 . 

15. 17.34.If ,riL£<T4044F1 , AS) .REVERT. 

15. 17. 34.  REVERT. 

15.17.34. REMAHE,TAPE2»T4044Ft. 

1 5 . 1 7 . 34 . ICO IN , OETf f N , , T4044F4 . 
15.17.3S.If.«0«.E8.m,LAB£LI. 

15. 17. 35.  BET .T4044F4/NA. 

15. 1 7.34. ELSE.LABEL1 . 

15. 17. 34.  ENDIF .LABEL) . 

15. 17.34.IF.F ILE<T4044F4,A$) .REVERT. 

15. 17. 34.  REVERT. 

15.17. 34 . RENAME , T APE3-T4044F4 . 

15. 1 7.34.6£TrMER6lTB/UN*BRYANTU. 

15.17.34. HER61TI. 

15.17.40.  CM  LUA-H  •  207111,  LOADER  USE!  344000 


15.17.41 .•••RET  EXTRACTION  USED 
15.17.43.  ••  INSERTIONS  BURINS  INPUT  ***mmmO 
15.17.43.  ••  OELETIONS  BURINS  INPUT  4MM44M0 
15.17.43.  ••  TOTAL  RECORBS  SORTED  •**•♦*•150 
15.17.43.  ••  INSERTIONS  BURINS  OUTPUT  •••*♦*•**0 
15.17.43.  ••  DELETIONS  BURINS  OUTPUT  •••*♦*•*•0 
15.17.43.  ••  TOTAL  RECORDS  OUTPUT  •*•*•<•158 

15.17.43.  ••  MERSE  ORBER  USED  •*«•*•*  11 

15.17.43.  4*CNB  SORT  RUN 
15.17.43.  END  MER6IT 


15.17.43.  050200  MAXIMUM  EXECUTION  FL. 

15.17.43.  O.Off  CP  SECONDS  EXECUTION  TINE. 

13.17.43. RENANE,T4044»TAP£4. 
1S.17.44.SAVEBI.T4044. 

15. 17. 45. REPLACE, T4044. 

15. 17. 44.  REVERT.  SAVE  < INDIRECT i 

15.17. 44 . CHANSE , T4044/CT-PM . 

IS.  1IL.44.BATf  ILE. 

•REVERT.  CCL  . 


gure  2.10:  An  Ix**plt  o t  tht  Intsractlvt  of  ths  MBRGIT  Proctdurt 
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continuation  of  previous  Fig 


Input  to  the  STATS  procedure  consists  of  a  title  card: 

Card  Columns  Description 

1  1-80  Title 

An  example  of  a  deck  setup  is  given  in  Figure  2.12.  This 
deck  is  available  as  CCSDS3/UN=BRYANTW. 


2.5.6  SDS  Time  Series  Plots 

It  is  possible  to  obtain  pen  plots  of  the  time  series  at  any 
point  in  the  normal  SDS  processing  cycle  (see  Reference  2) 
using  the  SDSPLOT  procedure.  Time  series  plots  can  be 
obtained  after  processing  by  DISC,  SAN,  FIL,  or  FFX.  The 
default  is  FIL.  Figure  2.13  shows  a  typical  time  series  plot 
after  processing  by  FIL.  Input  to  the  SDSPLOT  procedure 
consists  of  the  normal  SDS  processing  input  (Section  2.5.3) 
plus  two  other  parameters,  the  number  of  seconds  per  inch  and 
the  maximum  absolute  value  for  phase. 

Input  to  SDSPLOT  procedure: 

Record  1  -  Normal  SDS  input  (see  Section  2.5.3) 

Record  2  -  Card  Description 

1  SCPRIN,  PMAX 

where  SCPRIN  is  seconds  per  inch,  and  PMAX  is  maximum 
absolute  value  for  phase.  Both  quantities  are  REAL  and  in 
free  format  (i.e.  list  directed).  An  example  of  a  deck 
setup  for  SDS  time  series  plots  is  given  in  Figure  2.14. 

This  deck  is  available  as  CCSDS4/UN=BRYANTW. 


I 


•'.T(JP  Of  fu. o 
BRYANT .  SOS  STATISTICS 
USER . PRYANTU, PRYANTU . 
CHARGE, 3501 .4643. 

SET  ,f’ROCriL/UH»PRYAHTU . 
PESIH. STATS., T4044.SUA4044. 
EXIT. 

COR 

TRQHSO  TEST  CASE 
<B0TT0«  OR  fUE> 


Figure  2.12  Example  of  Dock  Setup  for  End  of  the  Month  Statistics 


2.5.7  Dumping  SDS  Field  Tapes 


It  is  possible  to  examine  the  data  on  SDS  field  tapes  using 
the  SDSDUMP  procedure.  The  data  are  converted  to  intensity 
and  phase,  but  no  detrending  or  correction  for  receiver 
artifacts  has  been  applied.  Input  to  the  SDSDUMP  procedure 
consists  of  the  normal  SDS  input  (See  Section  2.5.3).  An 
example  of  a  deck  setup  for  SDSDUMP  is  given  in  Figure  2.15. 
This  deck  is  available  as  CCSDS5/UN=BRYANTW.  Bear  in  mind 
when  dumping  files  that  each  8  second  time  block  produces  2 
pages  of  output. 


•<T0P  OF  f  ILE) 

BRYANT.  SDS  TINE  SERIES  PLOTS 
USER, BRYAN TU.BRYANTU. 

CHAR6C ,5504 .6670. 

SET  ,PROCnL/UN*BRYANTU. 

SET  ,TAPE19*T194<M4/UM*BRYANTU. 

LABEL . TAPE7.VSN*H91 ,NT,B*PE,F*S,LB«KU,PO*R. 
BE6 IN , 50SPLGT , ,FIL. 

E*IT. 

EOR 

TK4044 
1 ,4.1.84 
61 .21 ,31.0 
61,21 .43.0 
EOR 

60 .  .*>0 . 

<P*0T  TON  or  FILE) 


Figure  2.14  An  Bxsaple  of  Deok  Setup  for  SDS  Tine  Series  Plots 
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3.0  Study  of  Magnetic  Fields 


AFGL’s  research  emphasizes  investigations  of  the  ionosphere 
and  magnetosphere.  In  this  environment  Earth's  magnetic 
field  is  recognized  as  a  key  controller  in  the  global 
distribution  and  movement  of  particles,  including  pitch  angle 
dependencies,  mirroring,  drift,  and  dissipation.  At  AFGL, 
research  is  sponsored  on  satellites  such  as  DMSP,  SCATHA  and 
CRRES,  where  experimental  instruments  measure  the  particle 
environment  and  plasma  interactions.  This  chapter  describes 
the  installation  and  investigation  of  magnetic  field  models 
for  support  of  related  AFGL  research. 

3 . 1  Magnetic  Field  Models 


Below  a  few  earth  radii.  Earth's  internal  magnetic  field  is 
the  predominant  component.  The  internal  field  is  represented 
through  the  magnetic  field  potential  V  which,  at  any  point  in 
space,  is  expressed  in  a  spherical  harmonic  expansion: 


N  n 

V  =  a  2  (a)"*1  2  (g"  cos  m«l»  +  h"  sin  m<t>)  P™(sin  x) 

n=l  V  r /  m=0 


where 


and 


a  =  radius  of  earth,  $>  *  longitude,  x  =  latitude 


P„(sin  x)  =  associated  Legendre  polynomial 

The  vector  components  of  the  field  at  the  point  of  interest 
are  given  by: 


Br  =  -6V  =  -Z  (the  vertical  field  component) 

sr 

Bx  =  -Bq  =  -1  5V  =  +X  (the  North  field  component) 
r  sx 

Bo  =  1  SV  =  +Y  (the  East  field  component) 

r  sinx  s<t> 


The  external  magnetic  field  modifies  the  internal  field 
further  from  Earth.  Since  this  contribution  is  due  to  ring, 
tail,  and  magnetopause  currents,  it  must  be  expressed  in 
solar  oriented  coordinates. 

The  Earth-centered  dipole  magnetic  field,  based  on  a  1/r2 
potential,  is  the  simplest  but  still  useful  representation. 
The  Mead(1)  model  adds  a  symmetric  compression  as  well  as  a 
single  azimuthal ly  asymmetric  term  to  the  dipole  field,  and 
the  Mead-Williams<2)  model  includes  a  term  for  tail  currents; 
they  are  suitable  for  qualitative  studies  of  the  outer  zone. 
Full  scale  models  may  be  static,  ie.  with  no  model  for  geo¬ 
magnetic  variations,  such  as  the  Olson-Pf itzer  BFLD(3)  model 
used  for  SCATHA  studies;  or  they  may  allow  variations  due  to 
geomagnetic  activity,  eg.  the  Tsyganenko-Usmanov(4)  model. 


3.2  Magnetic  Field  Package  MFP 

More  than  10  years  ago,  AFGL/LCY  developed  and  integrated  a 
geomagnetic  field  package'5’  MFP,  for  general  use  at  the  Lab. 
The  routines  provide  B  and  L  values,  and  are  capable  of  field 
line  tracing.  A  spherical  harmonic  field  model  such  as  GSFC 
or  IGRF  is  input  as  Gauss  normalized  coefficients  repres¬ 
enting  the  zonal  and  tesseral  harmonics.  This  package  has 
served  most  of  the  lower  altitude  research  needs,  where  the 
external  field  is  negligible.  It  has  the  capability  of 
calculating,  for  a  given  geodetic  or  geocentric  position: 

o  Field  components  including  total  field 
o  Field  at  equator  crossing  &  conjugate  point 
o  L-shell  value 

o  Latitude  &  longitude  of  conjugate  point 

Coordinates  &  field  values  at  North  &  South  crossings 
of  a  specified  altitude 


o 
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A  separate  routine  LINTRA  performs  the  field  line  tracing  to 
specified  altitudes,  using  routine  FIELD  from  FIELDG,  but 
is  not  used  for  L-shell  calculations. 


3 . 3  Magnetic  Field  Package  BFLD 

The  magnetic  field  code  BFLD  was  developed  by  Aerospace  Corp. 
primarily  for  the  SCATHA  satellite  orbit  environment.  It 
contains  the  Barraclough<3)  model,  epoch  1975,  for  the  main 
internal  field,  and  an  Olson-Pf itzer  model  for  the  field  due 
to  external  currents.  The  external  field  is  calculated  using 
coefficients  that  are  functions  of  the  dipole  tilt,  and  this 
routine  BXYZMU  is  called  only  beyond  2  Re.  For  determining 
L-shell,  field  lines  are  traced  to  the  conjugate  point  using 
a  Runge  Kutta  4th  order  method,  and  the  invariant  integral  I 
is  computed  using  Simpson's  rule.  L  is  then  empirically 
computed  from  I  and  the  magnetic  field  strength. 


\ 

■ 

is 


3 . 4  Magnetic  Field  Package  Tsyqanenko-Usmanov  Model 
Coordinate  Systems  and  Rotations 

The  T-U  model  adds  Kp  dependence  to  the  Ring,  Tail,  and 
Magnetopause  contributions  of  the  external  field.  These 
components  are  formulated  in  different  coordinate  systems, 
all  Earth  centered.  Therefore,  the  program  initially 
calculates  the  rotation  matrix  transformations  that  will  be 
applicable.  Five  basic  coordinate  systems  and  two  aberrated 
(optional)  systems  are  involved: 


1 


Tsyganenko-Usmanov  Model  ( cont '  d ) 


Coordinate  System 

Primary  Axis 

Secondary  Axis 

Geographic  (G) 

x:  Greenwich, 

z:  North 

Uniform, 

free  rotation 

precession 

assumed 

Dipole  (D) 

zd:  N  mag.  pole. 

xd:  Geog.  meridian  of 

N.  Mag.  pole 

Celestial  (C) 

zc:  North  (=z). 

xc:  1st  point  of  Aries 

Solar  Ecliptic  (S) 

zs:  23.4333°  incl. 

xs:  points  to  sun 

to  ecliptic 

(Keplerian  orbit) 

Solar  Magnetospheric 

xsm:  towards  sun. 

zsra:  North  in  plane 

(SM) 

of  xsm,  zd 

Note:  zsm  rocks  111.2  daily  about  xs,  and  can 

be  inclined  as  much  as  34.6°  (23.4+11.2) 
to  the  ( xa ,  zg )  plane  at  equinox. 


OPTIONAL  (aberrated  coordinates) 

Solar  Wind  (SW)  resembles  Solar  Ecliptic  (S)  coordinates, 

but  xsw:  points  to  solar  wind;  e.g.  400  km/s 
solar  wind  and  20  km/s  Earth  orbital 
velocity  gives  an  aberration  angle  of  4.3°. 


Solar  Wind  Mag.  (SWM)  similarly  corresponds  to  SM. 


The  magnetic  vector  potential  is  represented  by 

A<t>  =  4  B03  p02  p2  [  z  +  p  +  4p0  ]'3/2 

in  the  cylindrical  geodipole  system  (p,«t»,z)  with 
2pc  =  the  loop's  radius.  B0  and  p0  are  Kp  dependent  model 
parameters  B0  =  -12.5  nT  (lowest  Kp)  to  -48.1  (highest  Kp) 
but  pQ  is  held  constant  at  4Re. 

In  terms  of  normalized  variables  p'  =  p/pQ  and  z'  =  z/pc 
this  gives 

Bz  =  4B0  ( 2z ' 2  +  8  -  P,2)/u 
Bp  =  2B0  p'z'/u 

2  2  5/2 

in  Dipole  coordinates,  where  u  =  (p'  +  z'  +4) 

Note  that  this  field  vanishes  at  z'  =  0,  p'  =  23/2. 


Tail  Current  Magnetic  Field 


An  infinite  sheet  of  current  filaments  is  assumed  parallel  to 
the  y-axis,  "hinged"  or  offset  a  distance  zD  on  the  z-axis 
due  to  the  dipole  tilt.  Kp  dependent  model  parameters  define 
the  current  per  unit  length  of  the  sheet  in  a  linear 
relationship: 


I(xQ)  =  c/2ir  (Bn  +  ABxo  -  x_) 

S 

where  S  =  xN-xF  and  N  and  F  refer  to  the  near  and  far  limits 
of  integration.  The  far  limit  is  fixed  at  50  Re. 

A  shaping  function  multiplier  [1  +  (y/Ay)2]"1  is  included  to 
make  the  current  flow  lines  curve  towards  the  sunward  edge  and 
approximate  the  ring  current  pattern.  A  closed  form  solution 
for  Bx,  By  and  Bz  has  been  derived  in  SM  coordinates. 


Magnetopause  Current  Magnetic  Field 


Bx ,  By  and  Bz  are  expressed  in  SM  coordinates  as  polynomial 
functions  of  y,  z,  exp(x/Ax),  and  sinY.  Ax  ranges  between  15 
and  20  Re  and  characterizes  the  gradient  of  the  external 
magnetic  field;  is  the  tilt  angle.  18  Kp  dependent 
parameters  are  the  coefficients  in  these  functions,  but  4  are 
mutually  related  due  to  the  imposition  of  the  Div  B  =  0 
condition. 


Organization  of  the  GTUBFD  program 


A  number  of  control  parameters  are  entered  in  an  array,  and 
activate  various  options  (or  defaults)  throughout  the 
program.  A  flow  chart  shows  the  steps  and  routines  employed 
when  x,  y,  and  z  are  input  in  SM  coordinates. 


Input 
x,y,  z 

( SM  coordinates) 


Routine  BFELD 


one-time  call 
unless  reactivated 


READY  (fixed  values) 
CONSTU  (Kp) 


TRANS3  ( rotations ) 


Routine  POINT 


Rotation 


ROTATE  coordinates  SM-G 
NSCALE  Truncate  Harm.  Expan. 
MAINB  Main  field 
ROTATE  B  G-SM 


ROTATE  coordinates  SM-D 
RING  B* 

ROTATE  Bj,  D-SM 
Add  Bj,  to  B 


Note: 

If  SWM  coordinates 
were  desired. 
Rotations  20,  19,  21, 
and  22  would  have 
been  used  instead 


TAIL  Bt  in  SM  coordinates 
MPAUSE  Bp,  in  SM  coordinates 
Add  Bx  and  B*  to  B 


Output 
Bx,  By,  Bz 
in  SM  coordinates, 
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4.0  Radiation  Belt  Studies 


The  high  energy  charged  particles  in  the  earth's  radiation 
belts  threaten  successful  operation  of  microelectronic 
circuits  of  satellites  operating  in  that  region.  AFGL  has 
conducted  research  on  the  radiation  belts  with  the  DMSP(1) 
satellites  at  low  altitudes  and  the  SCATHA<2)  satellite  at 
geosynchronous  altitudes.  In  this  section  we  summarize  the 
software  developed  for  analysis  of  SCATHA  data  and  display  of 
radial,  pitch  angle,  and  energy  distributions  of  trapped 
particles,  and  display  of  measured/model  magnetic  field 
comparisons. 

The  SCATHA  SC-5  instrument  measured  electron  and  ion  fluxes 
in  the  energy  ranges  50  eV  to  1  MeV  and  50  eV  to  7  MeV, 
respectively  (2) .  The  SC11  instrument  measured  the  magnetic 
field.  The  satellite  orbited  in  the  L-shell  range  5.3-8. 5  RE 
in  a  near-geosynchronous  orbit.  Software  was  developed  to 
form  and  display  orbital  data  bases  consisting  of  averaged 
detector  counts  by  L-shell/pitch  angle  bins  and  averaged 
magnetic  field  values  by  L-shell  bins,  with  the  corresponding 
ephemeris  parameters  appropriate  to  each  L-shell  bin. 
Subsequently  these  data  bases  have  been  used  in  a  study(3)  of 
adiabatic  variations  of  the  trapped  particles  for  the  purpose 
of  developing  schemes  which  might  be  used  in  the  analysis  of 
CRRES(4)  data. 

4 . 1  Statistical  Analysis 

Six  electron  energy  channels,  from  23  keV  to  1  MeV,  and  11 
ion  channels,  from  23  keV  to  5  MeV,  were  chosen  for  this 
analysis.  The  data  for  these  channels  are  available  at 
one  second  intervals  as  counts  per  .2  second,  along  with  the 
pitch  angle  at  6°  resolution  (satellite  rotates  at  360 
deg/min).  At  the  beginning  of  the  data  tape  (one  tape  per 
day)  are  header  records  indicating  the  detector  energies. 


geometric  factors,  and  degradation  factors.  These  are 
defined  so  that  the  fluxes  can  be  computed  from  the  counts  by 

F  =  GDC 

where  F  is  the  flux  in  particles/( cm2-sec-sr-eV ) ,  G  is  the 
geometric  factor,  D  is  the  degradation  factor,  and  C  is  the 
counts  per  .2  second.  The  degradation  factors  account  for 
variation  (degradation)  in  detector  sensitivity  in  time,  and 
effect  only  the  ESA  channels  below  100  keV.  Periodic  updates 
of  the  degradation  factors  occur  on  the  data  tape  as  short 
records  (i.e.  they  may  be  detected  in  processing  by  their 
short  length).  An  ephemeris  file  for  each  orbit  accompanies 
the  data.  This  includes  magnetic  field  parameters  based  on 
the  Olson-Pf itzer  quiet  model. <5) 


Program  SC5STAT  was  developed  to  process  this  data.  Data  are 
sorted  into  L-shell  bins  .05  RE  wide.  Each  L  bin  is  in 
general  encountered  twice,  in  the  ascending  and  descending 
half  of  the  orbit,  respectively.  These  are  treated  as 
separate  bins,  and  so  labeled.  The  data  for  a  particular 
such  bin  will  thus  be  encountered  in  one  continuous  period  in 
time.  Thus  the  data  can  be  processed  one  bin  at  a  time. 

Using  the  ephemeris  file,  SC5STAT  detects  which  L  bin  the 
satellite  is  in  initially  and  begins  accumulating  the 
necessary  statistical  sums  for  that  bin.  When  the  satellite 
leaves  that  bin,  these  sums  are  stored  on  a  temporary  file, 
and  the  process  repeated  for  the  next  L  bin.  Within  each  L 
bin,  the  data  are  sorted  into  the  10  deg.  pitch  angle  bins, 
and  separate  sums  accumulated  for  each,  along  with  the  count 
of  the  number  of  samples.  When  the  end  of  the  data  is 
reached,  SC5STAT  executes  a  second  pass,  in  which  the 
statistical  sums  are  converted  to  averages  and  standard 
deviations.  The  relevant  ephemeris  values  for  the  L  bin  are 
computed  by  interpolation  of  the  data  on  the  ephemeris  file 
at  the  average  universal  time  of  the  bin.  The  output  record 
for  each  L  bin  then  consists  of: 


rsr 


Average  universal  time; 

Total  number  of  readouts; 

L-shell  value  at  center  of  bin; 

Ascending/descending  flag  ( +1  if  ascending, 

-1  if  descending); 

Altitude  (km); 

Magnetic  latitude,  long.,  local  time; 

Latitude  and  local  time  in  SM  and  GSM 
coordinate  systems; 

Model  magnetic  field  vector  (SM)  and  magnitude; 

Magnetic  field  magnitude  at  equatorial  crossing 
of  field  line; 

Solar  position  vector  (ECI  coordinates); 

Dipole  moment  (ECI  coordinates); 

Mean  counts,  standard  deviation,  and  number  of  samples 
for  each  detector  for  each  pitch  angle  bin. 

The  calibration  factors  (geometric  and  degradation  factors) 
described  above  are  copied  to  a  separate  file  for  subsequent 
conversion  of  the  data  from  counts  to  fluxes. 

Program  MAGBIN  similarly  processes  the  magnetometer  data, 
producing  averages  for  the  same  L  bins  as  for  the  particle 
data  just  described.  For  each  L-bin  the  output  record 
consists  of: 

Average  universal  time; 

L-shell  value  at  bin  center; 

Mean  measured  and  model  magnetic  field  vector  components; 
Measured  and  model  magnetic  field  magnitudes; 

Measured  and  model  magnetic  pressure; 

Standard  deviations  of  measured  vector  components; 

Program  MAGMRG  combines  these  results  with  the  particle  data 
base  to  produce  a  single  particle  -  magnetic  field  data  base. 


Several  programs  were  written  for  display  of  the  data. 

Program  PLOTFIL  extracts  the  40°  and  90°  pitch  angle  data  for 
an  ascending  or  descending  portion  of  an  orbit  and  creates  a 
separate  file  of  this  data  suitable  for  plotting  with  program 
SUATEK(6) .  The  data  bases  for  two  consecutive  days  must  be 
used  since  the  requested  ascending  or  descending  portion  may 
not  lie  within  one  day.  If  the  data  requested  extends  past  a 
day  boundary,  the  bin  at  the  end  of  the  first  day  may  be  for 
the  same  L  as  the  first  bin  for  the  second  day,  i.e.  that  L 
bin  is  continued  into  the  second  day.  In  this  case,  the 
means  of  the  combined  bin  are  computed  as  the  weighted  means 
of  the  separate  portions,  using  the  sample  counts  that  had 
been  saved  for  the  separate  portions. 

Program  SIPLOT  is  a  specialized  version  of  SUATEK  developed 
to  stack  plots  sideways  for  up  to  11  detectors.  (Figures  4.1 
and  4.2)  For  each  detector  the  40°  and  90°  counts  are  placed 
in  a  single  panel.  A  separate  version  of  SIPLOT  similarly 
displays  the  measured  and  model  magnetic  field  vector 
components,  magnetic  field  magnitudes  and  magnetic  pressures. 

Program  SPECFIL  develops  spectral  plot  files  for  a  specified 
L  bin,  or  specified  group  of  selected  L  bins,  for  a 
specified  detector  and  pitch  angle.  The  data  on  this  file 
can  be  plotted  with  SUATEK.  Either  the  flux  or  the 
distribution  function  may  be  plotted.  The  flux  is 
computed  as  described  previously.  If  the  distribution 
function  is  desired,  it  is  computed  from: 

D  =  FXQ(E)/E, 

where : 

F  =  flux; 
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-1617  for  electrons,  5.45x10s  for  ions; 


E  =  energy; 

Q(E)  =  1/ 1 1+E/ ( 2E0 ) ] ; 

E0  =  rest  mass  energy,  5.11x10s  eV  for  electrons, 
9.4x10®  eV  for  ions; 

In  the  above,  the  flux  is  in  particles/ ( cm2-sec-sr-eV ) ,  the 
energies  are  in  eV,  and  the  distribution  function  is  in 
sec3/km6.  The  factor  Q(E)  is  a  relativistic  factor  which 
differs  significantly  from  unity,  for  the  energies 
considered,  only  for  electrons.  Otherwise  the  distribution 
function  given  here  is  identical  to  that  given  for  velocity 
space  in  Ref.  2.  Since,  due  to  relativistic  effects,  the 
momentum  distribution  function  is  not  related  to  the  velocity 
distribution  function  by  simple  multiplication  by  a  constant, 
we  have  chosen  to  use  the  relativistic  momentum  distribution 
function,  multiplied  by  the  cube  of  the  rest  mass,  so  that  it 
is  in  the  same  units  as  the  velocity  distribution  function 
quoted  in  Ref.  2,  and  reduces  to  it  when  relativistic  effects 
are  negligible. 

Spectral  analysis  programs  have  been  developed  in  Fortran  IV 
on  the  CYBER  to  display  measurements  from  SCATHA  in  color  on 
the  Tektronix  4100  series  terminals.  Particle  measurement 
and  ephemeris  data  are  read  from  the  SC5STAT  output  data 
tape,  interpolation  is  performed  and  the  data  are  arranged 
into  files  appropriate  for  plotting  using  color  graphics 
programs  to  interface  with  the  Tektronix  terminals.  The 
spacecraft  measurements  corresponding  to  electron  and  ion 
flux  distributions  in  space  are  thus  displayed  in  a  form 
ready  for  photographing. 

Figure  4.3  shows  a  sample  electron  spectrogram  generated  by 
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this  software.  In  this  case  the  particle  measurements  used 
cover  seventeen  energy  channels  and  ranges  from  5.5  to  8.5 
Re.  The  data  are  divided  into  60  L-bins,  each  equal  to  .05 
Re.  Each  L-bin  is  further  divided  into  10°  pitch  angle  bins. 
The  flux  value  for  the  particles  in  each  energy  channel, 

L-bin  and  pitch  angle  bin  is  read  from  the  SC5STAT  output 
data  base  and  is  stored  in  a  direct  access  file.  This  data 
are  then  arranged  for  graphics.  Each  pixel  in  the  plot  must 
be  assigned  a  color  to  display  the  flux  at  that  point.  The 
y-axis  on  each  plot  represents  the  particle  energies,  and  the 
x-axis  provides  sequential  particle  binning  data,  i.e.  all 
pitch  angle  bins  per  L-shell  bin.  Since  there  are  only  six 
electron  energy  channels  and  11  ion  energy  channels,  and  the 
resolution  of  the  energy  axis  is  90  pixels,  the  flux  values 
between  these  channels  are  interpolated  to  obtain  a  flux 
value  for  each  of  the  90  pixels.  The  entire  data  set  is 
displayed  as  two  sets  of  three  plots  --  one  for  electron 
measurements  (Figure  4.3)  and  the  other  for  ion  measurements. 
Each  set  of  three  plots  range  from  5.5  to  6.5  Re,  6.5  to  7.5 
Re,  and  7.5  to  8.5,  respectively.  The  x-axis  for  each  plot 
is  380  pixels  long,  so  there  is  a  pixel  for  each  flux  value 
in  each  L-bin  and  pitch  angle  bin,  therefore  no  interpolation 
between  flux  values  is  done  in  this  direction.  After  a  flux 
value  for  each  pixel  is  calculated,  15  flux  bins  are  defined 
over  the  entire  range  of  fluxes  and  a  bin  number  replaces 
each  flux  value.  The  bin  numbers  represent  the  color  index 
number  which  is  associated  with  each  pixel.  This  color 
index,  ranging  from  0  to  15,  defines  the  color  which 
represents  each  flux  bin.  (The  color  index  for  black  is 
zero.  This  index  represents  when  there  is  no  data  or  when 
the  flux  value  equals  zero. )  The  final  data  set  of  flux 
color  index  numbers  are  then  stored  in  a  direct  access  data 
files  to  be  used  with  the  plotting  program,  SPECPLT. 

Program  PTCHFIL  performs  a  function  similar  to  SPECFIL,  with 
respect  to  pitch  angle  distributions,  i.e.,  it  prepares  plot 
files  for  the  generation  of  pitch  angle  distribution  plots. 
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5. 


Investigation  of  Auroral  Electron  Precipitation  Data 

Electron  precipitation  flux  measurements  are  made  on  board 
the  DMSP  series  of  polar  orbiting  satellites  that  are  flown 
at  around  840  km  altitude.  The  data  bases  are  preprocessed 
and  made  available  for  a  number  of  research  investigations. 
Radex  personnel  have  conducted  computational  and  analytical 
work  for  statistically  interpreting  the  data,  and  providing 
estimates  for  observational  planning  purposes  of  instruments 
that  might  be  placed  on  shuttle  or  other  spacecraft  missions. 
For  these  studies  it  was  determined  by  the  researchers  that 
the  midnight  eguatorward  auroral  boundary,  which  can  be 
calculated  in  real  time  from  the  satellite,  provided  a 
superior  measure  of  geomagnetic  activity  which  should  be 
used  for  correlation. 

A  different  line  of  endeavor  consisted  in  developing  analy¬ 
tical  models  of  four  auroral  properties  viz.  electron  energy 
and  number  flux,  and  the  Peterson  and  Hall  conductivities. 
This  effort  led  to  functional  fits  using  spherical  harmonic 
and  Epstein  functions  requiring  a  relatively  small  number  of 
coefficients.  Through  contour  plots,  the  results  were  shown 
to  match  the  data  extremely  well,  and  could  be  used  both  for 
editing  irregularities  in  the  data  base,  and  for  prediction. 

These  investigations  have  required  accessing  large  data 
bases,  and  progressively  accumulating  and  managing  the 
reduced  data  for  the  analysis  and  presentation  purposes. 

Both  the  probabilistic  and  the  analytical  modeling  efforts 
were  successfully  completed,  using  F2,  F4,  F6,  F7  DMSP  data 
bases  on  individual  studies,  and  additional  runs  were  then 
requested  and  completed  on  other  data  bases. 


5 . 1  Probabilistic  Modeling  of  Auroral  Electron  Precipitation 
This  work  was  initiated  with  the  goal  of  improving  the 


quality  of  the  observation  of  bright  discrete  auroral  events 
through  the  calculation  of  the  probability  of  occurrence  of 
such  events  as  a  function  of  magnetic  activity  and  of 
geomagnetic  location.  With  such  knowledge,  one  could 
concentrate  observational  activities  only  during  periods  of 
magnetic  activity  consistent  with  a  high  probability  of  sight 
and  could  focus  attention  during  these  periods  on  regions  of 
the  sky  with  the  highest  probabilities.  Two  large  scale 
studies  were  carried  out,  the  first  on  data  from  DMSP/F2  and 
DMSP/F4  over  a  period  of  about  300  days(1)  and  the  second  on 
data  from  DMSP/F6  and  DMSP/F7  for  the  years  1983  and  1984 
respectively . (2)  The  major  difference  between  these  two  studies 
was  the  way  in  which  magnetic  activity  was  introduced.  In 
the  first,  the  Kp  index  was  used  and  in  the  second,  the 
latitude  of  the  midnight  extrapolated  equatorward  boundary  of 
diffuse  aurora  was  used.  The  correspondence  between  Kp  and 
midnight  extrapolated  boundary  has  been  established. (3) 

This  second  study  came  about  because  it  was  realized  that, 
since  equatorward  boundary  is  a  quantity  calculated  in  real 
time  for  each  DMSP  polar  pass,<4)  it  would  be  more  reliable 
to  refer  to  it  when  making  decisions  as  to  when  to  look  for 
high  energy  events.  This,  coupled  with  the  fact  that  the 
second  study  drew  upon  almost  twice  as  much  data  and  was 
derived  from  higher  resolution  instrumentation,  means  that 
for  practical  purposes,  the  second  study  superseded  the 
first.  We  will  thus  concentrate  this  discussion  on  the 
second  study. 

The  data  is  derived  from  DMSP/F6  and  F7  satellites.  F6  was 
launched  in  December  1982  into  a  840  km  circular  polar  orbit, 
sun- synchronous  in  the  dawn-dusk  meridian  plane.  The  F7 
satellite  was  launched  in  November  1983  in  a  similar  orbit  in 
the  1030-2230  meridian.  The  orbital  periods  of  both  are 
approximately  101  minutes.  Although  both  are  sun-synchro¬ 
nous,  the  spatial  coverage  for  the  two  taken  together  is 
quite  large,  due  to  the  rotation  of  the  geomagnetic  pole 
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about  the  geographic  pole.  Figure  5.1  shows  the  coverage  for 
the  periods  used  in  this  study. 

The  SSJ/4  sensors  on  F6  and  F7  are  identical  in  design  and 
measure  the  flux  of  electrons  and  ions  in  20  energy  channels 
in  a  range  from  30  to  30,000  eV.  For  the  electrons,  with 
which  we  will  be  concerned  here,  this  is  accomplished  by 
means  of  two  cylindrical  plate  electrostatic  analyzers,  one 
for  low  energy  and  one  for  high.  These  return  a  complete 
electron  spectrum  once  per  second.  A  more  complete 
description  of  the  experiment  has  been  given<5)  along  with 
results  of  the  calibration  necessary  in  conversion  from 
counts  to  flux  level.  The  DMSP  satellites  are  non-spinning 
and  the  sensors  are  mounted  such  that  their  look  directions 
are  oriented  radially  outward  from  the  earth  at  all  times. 

The  data  base  is  maintained  on  the  AFGL  Cyber  system  as  a 
series  of  tapes.  With  one  exception,  each  is  a  multi-file 
containing  either  five  for  F7  or  six  for  F6  files.  Each  file 
itself  contains  one  half  month  of  data  from  the  first  to  the 
fifteenth  or  from  the  sixteenth  to  the  end  of  the  month. 

Each  block  on  a  file  contains  5000  60-bit  CDC  words,  except 
for  the  final  block  of  each  day.  This  is  a  variable  length 
block  with  a  zeroed  block  of  length  two  immediately  before 
it.  This  means  that  if  one  wishes  to  read  past  a  day,  one 
must  look  for  a  block  of  length  two  and  then  read  one 
additional  block  before  reading  the  first  block  of  the  next 
day.  Each  minute  of  data  comprises  an  entity  of  variable 
length  within  the  fixed  length  block.  A  set  of  ephemeris  is 
assigned  to  the  first  second  of  each  minute  interval.  This 
variation  is  due  to  the  fact  that  missing  data  within  a 
single  record  is  not  zero  filled.  The  arrangement  leads  to 
two  complications.  The  fact  that  ephemeris  for  latitude  and 
magnetic  local  time  appears  each  minute  means  that  it  must  be 
interpolated  by  some  means  more  accurate  than  linear  to 
obtain  the  precision  we  require.  This  was  done  by  the  law  of 
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Figure  5.1  Magnetic  Latitude  -  Local-time  Coverage  by 

DMSP/F6  and  F7  Satellites.  (Regions  not  covered 
are  shaded) 
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cosines  as  outlined  in  the  next  section. 

Also,  the  fact  that  each  minute  of  data  is  of  varying  length 
means  that  minutes  overlap  blocks,  making  the  logic  for 
decoding  complicated.  We  have  created  a  routine  which 
returns  a  fully  decoded  minute  as  well  as  the  ephemeris  for 
the  next  minute  each  time  called.  With  this  subroutine  in 
hand,  we  can  ignore  the  complexities  of  the  data  base  for  the 
remainder  of  the  work. 

Once  each  minute  of  data  is  sorted  out,  the  nine  bits 
representing  each  channel  and  the  bits  for  the  ephemeris  are 
separated  and  channel  counts  are  converted  to  flux  by  the 
calibration  factors  mentioned  above. 


5.1.1  Processing  Algorithms 


In  order  to  assign  a  midnight  extrapolated  boundary  to  each 
DMSP  pass,  it  was  necessary  first  to  carry  out  the  extrapo¬ 
lation.  The  boundaries  themselves  were  obtained  from  the 
Hardy/Holeman  boundary  code  in  a  series  of  files,  one  month 
per  file.<4)  These  were  extrapolated  to  midnight  local  time 
by  the  relation  developed  by  Gussenhoven  et.al.<6)  The 
equatorward  midnight  boundary  is  given  in  their  notation  by 


-  (  ^OM  -  <XM/wi^0i  )  +  °V|/ °*i^CGM 


where  A  is  the  intercept  and  «  is  the  slope  of  the  Kp 
dependence  fit  for  the  boundary.  OM  refers  to  MLT  0, 

0i  refers  to  the  measured  MLT,  and  Acai  is  the  measured 
boundary.  The  result  A„  is  the  measured  boundary  extrapo¬ 
lated  to  its  equivalent  latitude  as  if  the  MLT  had  been 
midnight.  The  slopes  and  Intercepts  are  given  in  Reference 
6. 


Each  boundary,  morning  and  evening  side,  was  calculated 


separately  in  the  Hardy/Holeman  code  and  thus  there  were  two 
for  each  pass.  We  could  use  only  one,  since  a  midnight 
boundary  had  to  be  assigned  to  an  entire  pass.  We  chose  to 
use  the  evening  boundary  since  this  is  expected  to  be  the 
sharpest  and  most  well  defined  of  the  two.  It  was  necessary 
to  edit  the  boundaries  to  remove  some  poorly  determined 
points.  This  process  consisted  of  three  tests.  The  first 
was  made  in  the  boundary  code  itself,  which  rejected  poorly 
determined  cases.  Secondly,  we  required  the  extrapolated 
midnight  boundary  to  agree  with  that  predicted  on  the  basis 
of  Kp  to  within  two  standard  deviations.  The  predicted 
midnight  boundary,  also  from  Reference  6,  is  given  by 

■  aom  +  <*mkP 

Lastly,  we  required  that  a  boundary  appear  as  part  of  a 
morning/evening  pair,  that  is,  that  both  the  morning  and 
evening  boundaries  of  each  pair  passed  the  above  tests. 

There  was  another  complication  that  arose  because  we  wished 
to  include  data  at  all  latitudes  above  50  degrees.  The  UT 
associated  with  the  boundaries,  however,  was  the  UT  at  the 
boundary  which  was  usually  well  above  50  degrees.  In  order 
to  determine  the  UT  when  DMSP  passed  50  degrees  on  the  way  to 
encountering  each  boundary,  we  used  the  following  algorithm. 
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Figure  5.2  Geometry  of  50*  Latitude  Crossing 


Referring  to  Figure  5.2,  let  A  and  B  be  the  two  known  points 
in  the  orbit  given  in  the  boundary  file,  i.e.  the  position 
of  the  morning  and  evening  boundary.  Let  a,  b  and  c  be  the 
internal  angles,  measured  to  the  origin  of  the  sphere.  A  and 
b  are  thus  the  colatitudes  of  the  known  points.  Further,  the 
angle  Y  is  given  by  the  difference  in  MLT's  of  the  two  known 
points.  The  law  of  cosines  gives  us 

cos(c)  =  cos(b)cos(a)  +  sin( b )sin( a )cos( Y ) 

We  also  know  that  the  angle  c  must  vary  linearly  in  time 
since  the  satellite  orbits  at  a  constant  rate  and  constant 
radial  distance.  In  order  to  solve  the  problem,  we  need  to 
find  the  c  for  which  a  is  50  degrees.  First,  we  calculate 
«  from  the  law  of  cosines.  The  value  of  «  will  always  be 
the  same  no  matter  how  large  c  becomes,  so  from  here,  we 
solved  iteratively  for  the  proper  c  to  give  a  value  of  50 
degrees.  We  could  then  find  the  crossing  time  from 

ut50  =  UTb  -  (UTB-UTA)C50/C 

A  similar  algorithm  was  used  to  interpolate  the  latitudes  and 
local  times  between  the  once  per  minute  ephemeris  given  in 
the  DMSP  tapes. 


5.1.2  Quantities  Calculated 

The  major  probabilities  of  interest  are  what  we  call  per  pass 
and  per  second  probabilities.  The  per  second  probabilities 
are  defined  as  simply  the  fraction  of  seconds  in  which  an 


event  in  a  given  energy  and  flux  bin  was  observed,  divided  by 
the  total  number  of  seconds  spent  in  a  particular  geomagnetic 
division.  The  per  pass  probabilities  are  a  little  less 
straightforward.  These  give  the  fraction  of  passes  through  a 
particular  geomagnetic  bin  in  which  a  particular  energy  flux 
level  was  achieved,  divided  by  the  total  number  of  passes 
through  the  geomagnetic  bin.  In  order  to  qualify  as  a  pass, 
it  was  decided  that  the  satellite  must  spend  at  least  five 
seconds  in  a  given  bin. 

The  processing  system  requires  several  passes,  as  diagramed 
in  Figure  5.3.  First,  the  program  MIDER  examines  the 
boundary  file,  editing  out  poor  boundaries,  and  creates  a 
second  file  of  boundaries  for  a  particular  half  month.  The 
program  KOUNTA  then  uses  this  file  and  the  DMSP  tapes  to 
produce  three  files.  The  first  is  an  accumulation  of  per 
second  data,  which  is  binned  over  fixed  geomagnetic  locations 
and  midnight  equatorward  boundaries.  Second,  all  records 
that  contain  counts  greater  than  30,000  are  excluded  from 
this  count  and  dumped  to  another  file  for  visual  examination 
and  editing.  Third,  a  file  is  created  which  is  a  record  of 
each  geomagnetic  bin  encountered.  This  file  contains  two 
words  for  each  bin  crossing  of  the  form 

DDDTTTTT, MMMLLTTXNN 

where  DDD  is  the  day,  TTTTT  is  the  UT  in  seconds,  MMM  is  the 
equatorward  midnight  boundary  times  10,  LL  is  the  latitude 
bin  from  one  to  twenty  and  TT  is  the  mlt  bin  from  one  to 
forty-eight.  X  is  the  maximum  flux  level  encountered  in  the 
bin  from  one  to  six  and  NN  is  the  total  number  of  seconds  in 
that  pass  through  the  bin. 

Moving  to  the  next  level  in  the  flow  chart,  the  program  DUMPR 
reads  and  edits  the  high  counts,  creating  the  file  High 


Counts  which  contains  those  deemed  good  data.  These  are 
injected  into  both  the  per  second  data  and  the  bin  crossing 
data  files  through  the  programs  called  KOUNTC  and  COMBINE 
respectively.  The  per  second  data  is  then  combined  to  yield 
the  final  per  second  probabilities.  Yet  another  program  is 
then  used  on  the  files  called  xnyyy,  the  records  of  bin 
crossings,  to  calculate  the  per  pass  probabilities.  The 
beauty  of  the  per  pass  probability  calculated  in  this  way  is 
that  one  can  partition  the  count  over  any  arrangement  of 
midnight  boundaries  or  over  any  arrangement  of  geomagnetic 
bins  without  recourse  to  the  DMSP  data  base  itself.  We  have 
used  this  system  for  both  these  variations  on  the  original 
theme . 

5.1.3  Results  and  Deliverables 

The  results  of  both  these  studies  were  of  numerous  forms. 

The  most  significant  were  contour  plots  of  probability. 

There  were  four  sets  in  all.  The  first  came  from  the  per 
second  probabilities,  calculated  over  midnight  boundaries 
chosen  to  correspond  to  the  thresholds  2,4  and  6  in  Kp. 

The  second  set  was  the  per  pass  version,  calculated  over  the 
same  boundary  divisions.  A  third  set  was  produced  calculated 
over  four  boundary  divisions  that  were  chosen  such  that  each 
division  contained  roughly  the  same  number  of  passes.  The 
fourth  set  was  calculated  over  five  divisions  selected  in  the 
same  way  for  more  resolution. 

In  making  these  plots,  it  was  necessary  to  fill  in  the  small 
regions  of  space  which  were  covered  either  poorly  or  not  at 
all  by  the  orbits  of  DMSP  F7  and  F6.  Additionally,  a  mild 
form  of  smoothing  was  performed  on  the  data  before  plotting. 
This  consisted  of  averaging  over  each  set  of  three  local 
times  for  a  constant  latitude.  A  typical  example  of  such  a 
contour  plot  is  shown  in  Figure  5.4. 
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Another  result  of  interest  was  histograms  of  the  duration  or 
spatial  extent  of  each  high  flux  event.  A  typical  example  is 
shown  in  Figure  5.5.  Tables  of  numerical  values  for 
probability  were  also  produced.  Finally,  the  secondary  data 
base  discussed  above  was  used  to  calculate  a  rather  different 
probability,  the  probability  of  sighting  an  arc  greater  than 
or  equal  to  a  given  intensity  anywhere  in  the  auroral  region 
during  a  given  pass  at  fixed  equatorward  boundary.  Seven 
divisions  of  boundary  were  selected  for  this  study  and  the 
results  are  shown  in  Table  5.1.  We  note  again  that  without 
the  secondary  data  base  that  this  work  produced,  such  a 
calculation  would  have  had  to  rely  on  the  DMSP  tapes,  which 
would  have  required  weeks  instead  of  a  few  hours  to  complete. 

A  by-product  of  these  two  studies  was  the  substantiation  of 
the  correlation  of  Kp  with  midnight  equatorward  boundary. 

This  correlation  was  borne  out  by  the  qualitative  similarity 
of  the  results  of  the  first  and  second  study. 

5.1.4  Further  Work 

Interest  has  been  expressed  in  the  initiation  of  other 
counting  studies  on  the  DMSP  data.  The  major  difference 
between  these  and  the  previous  ones  would  be  the 
concentration  on  the  overall  structure  of  the  aurora  in 
regard  to  energy  flux  and  average  energy  rather  than  to  focus 
only  on  high  energy  events  as  was  previously  done.  This 
could  be  carried  out  by  increasing  the  number  of  energy  and 
flux  bins  and  by  lowering  the  thresholds  so  that  the 
probability  became  distributed  over  several  bins  instead  of 
concentrated  in  the  lowest  of  the  bins.  Binning  by  average 
energy  and  flux  would  then  show  which  regions  of  the  sky 
exhibited  high  energy  electrons  as  well  as  a  high  level  of 
flux.  Results  of  such  a  study  would  be  of  interest  in 
understanding  the  detailed  structure  of  the  aurora  as  opposed 
to  determining  when  and  where  a  bright  event  might  occur. 
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Figure  5.5  Example  of  Flux  Event  Duration  Histograms 


This  effort  could  be  accomplished  relatively  easily  using  the 
existing  processing  system,  which  is  tailored  for  DMSP/F6  and 
F7  data,  since  the  problems  involving  access,  decoding  and 
assigning  the  equatorward  boundaries  to  each  polar  pass  have 
been  encountered  and  successfully  solved  in  the  study  above. 

5 . 2  Modeling  of  Averaged  Auroral  Properties 

Another  ?rea  of  attention  in  regard  to  the  investigation  of 
electron  precipitation  data  was  the  modeling  of  averaged 
auroral  properties  by  simple  functional  forms.  Four 
properties  were  studied;  the  integral  and  energy  flux  and  the 
Hall  and  Peterson  conductivities.  These  properties  were 
calculated  by  Hardy,  et.al.  from  DMSP  F2  and  F4  data.<7)  Their 
model  averaged  over  time  at  a  series  of  Kp  values,  using  48 
divisions  of  magnetic  local  time  and  30  divisions  of 
latitude.  A  total  of  15  months  of  data  was  used,  chosen  so 
as  to  provide  a  uniform  distribution  over  seasonal 
variations.  The  properties  to  be  fit  were  obtained  as  a 
matrix  of  1350  values  at  seven  integral  values  of  Kp  from 
zero  to  six.  Each  matrix  contained  data  at  48  magnetic  local 
times  and  30  geomagnetic  latitudes.  In  general,  the  data  was 
statistically  excellent,  although  occasional  irregularities 
did  exist.  One  purpose  of  this  modeling  was  to  smooth  out 
these  irregularities,  although  the  most  important  objective 
was  to  reduce  the  data  to  a  more  manageable  size  via  the 
functional  fits.  The  rationale  was  that  if  the  large  arrays 
could  be  reduced  to  a  small  number  of  coefficients,  the  data 


could  be  employed  on  microcomputers  in  the  field  for 
predictive  purposes. 

Overall,  the  properties  as  viewed  as  a  function  of  latitude 
at  fixed  local  time  resemble  two  straight  lines  that  converge 
smoothly  to  a  maximum  at  a  characteristic  mid-latitude.  This 
regular  pattern  was  a  considerable  aid  in  choosing  the  fit 
functions. 


5.2.1  Methods  of  Representing  the  Properties 


Our  first  attempt  at  the  creation  of  a  model  used  a  spherical 
harmonic  expansion.  A  software  system  was  created  to 
interpolate  the  values  onto  a  49  x  49  point  grid  and  to 
calculate  the  expansion  coefficients  coefficients  by  direct 
integration  of  the  product  of  the  function  to  be  modeled  and 
each  of  the  spherical  harmonics.  This  particular  grid  was 
chosen  because  it  gave  an  accurate  result  for  the  mean  square 
integral  of  all  spherical  harmonics  up  to  order  seven  and  yet 
was  not  too  time  consuming. 

The  representation  of  spherical  harmonics  used  was 

7  n 

l°g10f  (  T,  X  )  =  2  2  b™  Y£(T,X) 

n=l  m=-n 

where 

Yjj  =  cos(  2irT/24  )P„(  X  ) 

Y-m  =  sin(  2ttT/24  )P™(  X  ) 

and 

T  is  the  magnetic  local  time  in  hours 
X  is  (latitude  -  70)/20 

and 

is  the  calculated  expansion  coefficient 

b™  is  a  coefficient  that  makes  the  square  integral 
of  each  spherical  harmonic  unity. 


The  use  of  C  and  b  is  convenient  because,  when  normalized  by 
b,  the  expansion  coefficient  is  simply  the  value  of  the 
integral  of  the  spherical  harmonic  and  the  function  to  be 
fit.  Its  use  does  mean,  however,  that  one  must  carry  along 
all  the  b's  in  the  generating  routine.  To  avoid  this  in  the 
implementation  of  the  model,  we  have  reported  coefficients 
both  as  the  C  and  as  the  product  of  b  and  C  in  the  above 
equation. 

This  scheme  was  used  to  expand  the  base  ten  logarithm  of  the 
energy  and  number  flux  to  produce  a  total  of  64  coefficients 
at  each  value  of  Kp.  The  conductivities  were  not  so 
amenable  to  fitting  with  spherical  harmonics  and  this 
approach  was  discarded  in  their  case  in  favor  of  a  second  fit 
function,  Epstein  transition  functions. 

In  our  implementation,  the  Epstein  function  was  taken  to  be 

e(h)  =  r  +  Si( h-h0  )  +  (  s2-Si  )  lnCl-Si/s^xpC  h-h0 )  ]  /  [  l-s^Sj] 

where  h  is  the  latitude  and  e  is  comprised  of  four  parameters 
the  maximum  value  r,  the  latitude  where  the  maximum  occurs 
hQ,  and  two  slopes  sx  and  s2  on  either  side  of  the  maximum. 
Some  distance  from  hQ,  the  Epstein  function  is  nearly  a 
straight  line  with  slope  sx  or  s2.  Near  h0,  these  two  lines 
meet  in  an  arc  which  reaches  maximum  r  at  h0.  This  is  a 
special  adaptation  of  the  more  general  Epstein  function. <8) 

Fitting  was  done  in  two  passes.  First,  the  maximum  latitude 
and  value  were  chosen  from  the  data  and  a  least  squares 
optimum  value  for  sx  and  s2  was  calculated  for  all  points 
above  a  cutoff.  In  order  to  inspect  the  quality  of  these 
fits,  an  interactive  graphics  system  was  developed  to  page 
through  each  local  time  and  Kp,  overlaying  the  original 
data  and  fit  function.  If  a  fit  were  judged  visually  to  be 
something  less  than  the  best  possible,  any  of  the  four 
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parameters  could  be  modified  in  one  of  several  ways  and  new 
least  squares  fits  could  be  performed  over  variable  ranges  of 
latitude.  Such  modification  was  found  to  be  necessary  in 
about  twenty  percen1-  of  the  cases. 

When  fitting  was  complete,  we  were  left  with  192  coefficients 
for  each  level  of  Kp.  Since  the  local  time  dependence  of 
the  coefficients  showed  a  slow  and  regular  variation,  we  felt 
that  the  set  could  be  represented  by  expansion  of  the  Epstein 
coefficients  in  a  Fourier  series  in  local  time.  We  took  this 
series  to  order  six. 

6 

<x(T)  =  2  [  Cncos(  2irT/24 )  +  Snsin(  2irT/24 )  ] 
n=0 

where 

T  is  MLT  in  hours  and 
«  is  r  or  h0  or  sx  or  s2  above. 

5.2.2  Results  and  Deliverables 

The  main  result  of  this  work  were  the  sets  of  coefficients 
and  software  to  implement  them.  As  support,  contour  plots  of 
original  data  and  model,  and  histograms  and  plots  of  percent 
error  were  produced.  A  typical  contour  plot  is  shown  in 
Figure  5.6  and  a  histogram  of  fractional  deviation  is  shown 
in  Figure  5.7.  Also,  tables  of  coefficients  were  supplied. 

In  the  case  of  the  Epstein  functions,  these  coefficients  were 
doubly  useful  in  that  they  could  be  used  to  give  directly  the 
position  of  the  maximum  and  the  minimum  value.  This  study 
has  been  documented*9*  and  is  being  prepared  for  submission 
to  the  Journal  of  Geophysical  Research. (10) 
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Figure  5 . 6  Example  of  Contours  from  Original  Energy  Flux 
Data  (top).  Spherical  Harmonic  Model  (bottom) 
and  Fractional  Deviation  as  a  Function  of 
Value  (right) 
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6 . 0  The  AFGL  Interactive  Targeting  System 

6 . 1  Introduction 

The  AFGL  Interactive  Targeting  System  (AITS)  was  designed  to 
provide  pre-mission  and  on-orbit  planning  support  for  the 
CIRRIS-1A  project.  It  is  a  computer  based  system  which 
displays  shuttle  position  and  sensor  line  of  sight  in  various 
formats.  Color  graphic  displays  are  used  to  depict  the  world 
map,  the  celestial  sphere,  and  the  probability  of  auroral 
activity.  This  system  is  required  to  provide  researchers 
with  the  tools  needed  to  evaluate,  in  near  realtime,  the 
shuttle  and  the  sensor  orientation  during  data  gathering 
portions  of  the  orbit,  and  to  view  the  dynamics  of  the 
shuttle  attitude  and  the  auroral  oval  conditions.  It  also 
provides  information  on  sensor  line  of  sight  and  field  of 
view  for  earth  and  space  targets.  The  system  is  user 
friendly  and  provides  the  researcher  with  an  interactive 
menu  to  select  the  desired  computation  or  display. 

The  AITS  system  was  developed  to  support  the  CIRRIS-1A 
project,  however,  the  system  can  be  applied  toward  any 
orbiting  sensor  whose  pointing  direction  can  be  calculated. 
The  system  is  especially  useful  for  sensors  whose  pointing 
direction  can  be  controlled. 


6.2  General  Description 


6.2.1  System  Application 

The  upper  atmosphere  is  highly  dynamic  and  exhibits  infrared 
spectral  and  spatial  structure  which  can  interfere  with  Air 
Force  systems'  detection  of  infrared  targets.  The  purpose  of 
CIRRIS-1A  is  to  measure  and  characterize  the  spectral  and 
spatial  properties  of  the  earth  limb  atmosphere.  Specific 
objectives  include:  infrared  radiation  from  atmosphere 
airglow,  auroral  emissions,  targets  of  opportunity,  and 
effects  of  local  shuttle  contamination. 

CIRRIS-1A  will  measure  the  spectral  and  spatial 
characteristics  of  the  30-300  km  earthlimb  over  a  range  of 
latitudes,  day/night  conditions,  and  geomagnetic  activity. 
This  data  will  form  the  basis  of  an  earthlimb  model  of 
spectral  signatures  and  spatial  clutter  which  can  be  used  by 
designers  to  optimize  the  performance  of  their  operational 
systems.  Secondary  objectives  are  to  measure  targets  of 
opportunity  and  to  demonstrate  the  role  of  a  military  payload 
specialist  (Manned  Spaceflight  Engineer  -  MSE)  in  on-orbit 
shuttle  operations. 

To  obtain  this  data,  the  CIRRIS-1A  measurement  modes  include 
vertical,  horizontal,  and  stare  scans  of  the  earthlimb  using 
various  interferometer,  radiometer,  optical  filters,  and  scan 
pattern  combinations.  For  example,  an  earthlimb  scan  be 
entirely  pre-programmed  as  to  where  the  sensor  will  step 
through  a  sequence  of  tangent  heights  to  provide  a  vertical 
atmospheric  profile.  On  the  other  hand,  for  transitory 
observations,  the  payload  specialist  will  track  auroral  or 
other  targets  of  opportunity  using  a  low  light  level  TV 
( LLLTV )  co-aligned  with  the  sensor  telescope,  and  a  joystick 
controller  in  the  orbiter  aft  flight  deck. 


This  report  provides  a  functional  description  of  the  ground 
based  flight  support  computer  system.  The  intended  uses  of 
the  system  are  to  provide  the  following: 

1.  Provide  the  tools  required  for  the  AFGL  researcher  to 
effectively  manage  the  data  measurement  periods 
available  to  the  CIRRIS  payload  on  orbit. 

2.  Provide  support  to  the  payload  specialists  as  required 
for  the  targets  of  opportunity. 

3.  Provide  limited  prefilight  mission  planning. 

4.  Provide  limited  post  flight  ancillary  data  analysis. 


The  AITS  software  operates  on  a  Digital  Equipment  Corporation 
(DEC)  VAX  11-750  computer  system.  The  computer  system 
configuration  consists  of  a  456  Mega-byte  fixed  hard  disk  for 
database  storage;  a  9-track  magnetic  tape  drive  for  permanent 
archive  of  databases  generated  and  for  initial  loading  of  the 
AITS  and  system  software;  a  system  printer  for  hard  copy 
output  of  selected  computations;  color  graphic  terminals  for 
use  in  the  generation  of  mission  required  displays;  and  a 
color  ink  jet  hard  copy  unit  to  obtain  copies  of  the  graphic 
terminal  displays  for  archive  and  for  use  in  mission 
planning. 


6 . 3  Operations  Support  Functions 

The  AITS  software  system  provides  researchers  with  the 
capability  to  compute  and  display  information  relevant  to  the 
operation  of  probes  flown  on  the  space  shuttle.  Data 
displays  are  used  for  long  and  short  range  operations 
planning,  and  post  pass  evaluation  of  data  taking  conditions. 
The  processes  calculate,  for  the  CIRRIS-1A  payload,  the 
sensor's  line  of  sight  and  orientation  needed  to  view  target 
locations.  The  CIRRIS-1A  payload  is  a  multi-sensor  gimbal 
mounted  system  with  the  sensor's  orientation  requirements 
determined  by  the  data  gathering  operation  being  performed. 

A  description  of  the  processes  and  displays  and  their  use  in 
the  planning  function  is  given  below. 


6.3.1  Vehicle  Flight  Simulation 

This  procedure  provides  displays  of  shuttle  ground  tracks  on 
latitude  and  longitude  grids  with  continental  outlines.  This 
can  be  either  a  either  linear  or  polar  coordinate  display. 
Also,  the  Feldstein(1>  auroral  oval  and  Hardy/Gussenhoven*2-3  4) 
probability  contours,  described  in  Section  6.4.4  of  this 
report,  and  the  day/night  terminator  locations  at  0  km,  100 
km,  and  satellite  altitude  are  projected  onto  these  maps.  A 
GMT  and  MET  time  scale  is  shown  on  the  bottom  of  the  display 
along  with  a  one  day  mission  time  line  containing  planned 
sensor  modes  of  operation.  These  displays  can  be  used  by  the 
researchers  to  determine  ground  tracks  of  the  past  and 
upcoming  orbits;  present  time  and  position  of  the  vehicle; 
line  of  sight  fan  of  the  instrument  relative  to  the  auroral 
oval  and  probability  contours;  and  when  the  next  instrument 
operation  mode  will  occur. 


Figure  6.1  shows  the  linear  world  map  projection  display  that 
is  used  to  monitor  the  position  of  the  vehicle.  The 
day/night  terminators,  Feldstein  auroral  ovals, 
Hardy/Gussenhoven  contours,  and  the  vehicle's  ground  track 
are  recalculated  and  displayed  at  a  rate  of  every  5  minutes 
during  the  monitoring  process.  The  ground  track  position  is 
shown  on  the  terminal  display  in  white,  directional  arrows 
depict  where  the  vehicle  was  for  the  past  30  minutes.  The 
portion  in  red  with  directional  arrows  shows  where  it  will  be 
in  the  next  120  minutes.  The  point  representing  the 
vehicle’s  location,  at  the  displayed  time,  is  depicted  with  a 
flashing  image  of  the  vehicle  on  the  ground  track.  The  image 
of  the  vehicle  is  repositioned  every  minute  to  represent  the 
actual  movement  of  the  vehicle  as  a  function  of  time.  The 
shuttle  pointer  below  the  time  scale  shows  the  present  time 
and  the  instrument  operation  occurring. 

The  user  can  select  one  of  three  vehicle  attitude 
orientations;  Belly  to  Ram  ( BTR ) ,  Gravity  Gradient  South 
(GGS),  or  Gravity  Gradient  North  (GGN)  when  determining  the 
region  for  the  gimbaled  sensor's  field  of  view.  The  sensor's 
field  of  view  is  displayed  at  each  directional  arrow  as 
gimbal  scan  regions.  Each  gimbal  scan  region  is  displayed  on 
the  terminal  in  either  white  depicting  where  the  scan  region 
was  in  the  past  30  minutes,  or  red  depicting  where  the  scan 
region  will  be  in  the  next  120  minutes. 
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6.3.2  Ephemeris  Computations 


Orbit  parameter  computation  is  done  using  a  software  package 
derived  from  the  ephemeris  calculation  program  LOKANGL . <5-6) 

The  package  interpolates  between  orbital  elements  to  compute 
satellite  position,  or  extrapolates  from  an  orbital  element 
set  for  cases  of  thrust  or  lack  of  further  elements.  In 
interpolation  or  extrapolation,  the  package  uses  the  six 
elements  (a,  the  semi-major  axis;  e,  eccentricity;  N,  the 
right  ascension  of  the  ascending  node;  w,  the  argument  of 
perigee;  i,  the  inclination;  and  M,  the  mean  anomaly)  along 
with  their  derivatives  to  compute  sub-satellite  position. 

The  condensing  and  retailoring  of  the  LOKANGL  package  so  that 
it  is  included  in  the  software  system,  minimizes  the  need  for 
pregenerated  ephemeris  information  and  the  amount  of  computer 
storage.  The  major  advantage  of  these  independent  ephemeris 
computations  arises  when  online  updating  of  the  satellite 
position  elements  occurs.  Regeneration  of  data  files  to 
incorporate  the  change  is  unnecessary  since  the  routines 
access  the  orbital  position  elements  directly.  Initial 
orbital  elements  are  entered  prior  to  launch  and  are  updated 
as  they  become  available  during  on-orbit  operations. 


6. 3. 2.1  Trajectory 

Ephemeris  computations  performed  include  calculation  of 
vehicle  position  in  geodetic  and  geocentric  coordinates  and 
corrected  geomagnetic  coordinates.  Vehicle  position 
information  may  be  generated  in  the  form  of  hardcopy  listing 
or  as  a  graphical  display).  The  graphical  display  provides 
the  user  with  the  geographic  projection  of  the  position  of 
the  vehicle  for  the  time  span  selected.  Modified  mercator, 
linear,  or  polar  projections  of  the  world  map  may  be  used  to 
depict  the  geographic  outline  of  the  land  masses  and  the 
vehicle's  position. 
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6. 3. 2. 2  Station  Look  Angle 


The  station  look  angle  option  provides  the  user  with  the 
capability  to  compute  and  display  information  when  the 
shuttle  is  within  the  field  of  view  of  ground  stations. 

NASA,  SCF,  and  ADCOM  ground  station  coordinates  reside  in 
permanant  storage  and  may  be  selected,  or  the  coordinates  of 
any  ground  station  may  be  entered.  The  user  selects  the 
desired  minimum  elevation  for  the  station  viewing  of  the 
vehicle  during  the  time  span  selected. 

A  tabular  output  containing  start  and  end  times  of  viewing, 
and  the  maximum  elevation  time  at  each  station  is  generated 
(Figure  6.2).  Graphical  displays  of  modified  mercator, 
linear,  or  polar  projection  of  the  world  map  can  be  chosen  by 
the  user.  These  displays  contain  the  vehicle's  ground  track 
along  with  the  chosen  station's  conical  viewing  region 
determined  by  the  altitude  of  the  vehicle  and  the  minimum 
viewing  elevation  selected. 


6.3.3  Celestial  Support  Functions 


The  AITS  software  uses  various  displays,  symbols  and  colors 
to  represent  the  celestial  sphere  and  stellar  sources. 

Stellar  intensities  are  denoted  by  various  colors,  and  symbol 
sizes  within  a  color.  Each  display  shows  for  the  time  period 
selected,  the  celestial  ecliptic  plane,  and  the  planets.  The 
vehicle's  orbital  path  for  the  time  period  selected  is 
projected  onto  the  celestial  sphere  with  directional  arrows 
at  5  minute  intervals.  The  celestial  sphere  and  computed 
information  can  be  depicted  as  spherical,  aitoff,  or  linear 
display  projections.  These  displays  can  be  used  to  show  what 
the  instrument  will  view  in  the  celestial  background. 
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Figure  6.2  Station  Look  Angle  ( Example  Listing  #1) 


A  spherical  projection  of  the  celestial  sphere  consists  of 
approximately  one  fourth  of  the  sphere,  where  the  focal  point 
of  the  projection  is  user  determined.  The  central  region  of 
a  spherical  projection  supplies  the  user  with  an  actual 
representation  of  the  celestial  area  viewed.  This  type  of 
projection  allows  the  user  to  obtain  necessary  information 
about  the  LOS  or  FOV  of  the  sensor  as  it  tracks  along  a 
specific  region  of  the  celestial  sphere  during  that  orbit. 

Aitoff  and  linear  projections  allow  the  user  to  display  the 
complete  alestial  sphere  and  all  the  computed  information  on 
one  display.  An  aitoff  projection  minimizes  distortion  of 
the  polar  regions.  Both  types  of  displays  show  a  full  orbit 
projection  and  also  all  the  regions  that  the  sensor's  LOS  or 
FOV  encounter.  The  aitoff  and  linear  display  are  the  most 
useful  type  of  display  when  performing  gimbal  angle 
determination  for  a  celestial  target  sighting.  These 
projections  allow  the  user  to  view  all  the  gimbal  angles 
necessary  at  each  point  in  the  orbit. 
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6.3.4  Camera  Computations 


6.3. 4.1  Camera  View  Simulation  at  an  Instant 

The  AITS  system  process  for  Camera  View  Simulation  allows  the 
user  to  graphically  display  the  camera  field  of  view  ( FOV )  at 
a  given  instant  in  the  mission.  The  user  selects  the  time 
desired  along  with  the  vehicle  attitude  and  gimbal 
orientation  of  the  camera  to  be  used  for  the  simulation.  The 
graphic  display  contains  land  masses,  stars,  planets,  and  a 
simulation  of  the  auroral  boundary  that  would  be  encountered 
at  that  time  in  the  mission  (Figure  6.3).  The  FOV  is 
annotated  with  a  crosshair  that  marks  the  line  of  sight  (LOS) 
of  the  camera  and  the  angular  dimension  of  the  FOV.  The 
display  is  annotated  with  either  the  geographic  location  of 
the  LOS  tangency  or  the  intersection  with  the  earth.  In  the 
cases  where  the  camera  FOV  does  not  contain  any  portion  of 
the  earth,  the  display  is  anotated  with  the  celestial 
coordinates  for  the  LOS. 


6. 3. 4. 2  Camera  Celestial  Target  Tracking 

This  process  allows  the  user  to  select  a  celestial  target 
location  and  obtain  the  gimbal  angles  required  to  maintain 
the  target  within  the  camera's  FOV.  The  user  selects  the 
period  of  interest  in  the  mission  and  a  constant  vehicle 
attitude  for  that  period.  The  process  then  computes  the 
gimbal  angles  required  for  the  camera  to  view  the  target 
location  over  the  span  of  time  selected.  This  information  is 
shown  across  the  bottom  of  the  display  as  XY  graphs  (Figure 
6.4).  The  top  of  the  display  contains  the  celestial  region 
that  the  gimbal  tolerances  allow  the  camera  to  view. 


gure  6.3  Camera  View  Simulation  at  an  Ins 


The  actual  camera's  FOV  location  for  the  selected  time  is 
displayed  within  the  gimbal  region  with  a  red  outline.  The 
user  may  then  step  forward  in  time  to  view  the  camera's  FOV 
and  gimbal  region  at  another  instant  using  the  same  attitude 
and  target  parameters.  The  gimbal  region  display  contains 
the  stars  and  planets  viewed  and,  for  cases  of  earth 
intersection,  a  brown  spherical  representation  of  the  earth's 
region.  The  gimbal  angles,  vehicle  attitude,  and  the  time 
period  in  the  mission  are  stored  as  part  of  the  user  defaults 
during  this  process.  This  allows  the  user  the  option  of 
implementing  the  Camera  View  Simulation  process  to  generate  a 
camera  FOV  display  with  the  earth  land  masses,  stellar 
values,  and  auroral  region  encountered. 


6.3.5  Gimbal/Of f-Track  Calculations 

6.3. 5.1  Ground  &  Celestial  Target  Gimbal  Angle  Determination 

This  process  determines  the  gimbal  angles  required  to  view 
either  a  ground  target,  auroral  region,  or  a  celestial  target 
location  for  a  user  entered  constant  vehicle  attitude.  An 
auroral  region  is  defined  by  the  geometric  center  (centroid) 
of  a  selected  Hardy  Auroral  Probability  Contour.  The  results 
of  the  computation  performed  for  the  selected  span  of  time 
are  tabulated  on  print-outs  to  aid  the  researcher  in  the 
selection  of  gimbal  angles  to  use  with  other  processes.  The 
Local  Vertical  Local  Horizontal  ( LVLH )  angles  of  azimuth  and 
elevation  from  the  vehicle  to  the  target;  the  local  time  at 
the  target;  and  the  geodetic  position  of  the  vehicle  are 
tabulated  for  each  of  the  target  calculations.  LVLH 
elevation  is  referenced  positive  down  from  the  local 
horizontal  toward  the  earth.  LVLH  azimuth  is  positive  about 
the  local  vertical  from  the  negative  velocity  vector  on  the 
local  horizontal  plane. 


6. 3. 5. 2  Ground  &  Celestial  Target  LVLH  Angle  Determination 

This  process  determines  the  LVLH  angles  required  to  view  a 
ground  location  or  a  celestial  target.  The  user  supplies  the 
type  of  target  and  its  coordinates,  the  span  of  time  desired 
to  process,  and  the  type  of  output  product  desired.  The 
process  creates  either  a  tabulated  print-out  of  the  infor¬ 
mation  or  a  graphic  display  of  the  LVLH  angles  versus  time. 

The  print-out  values  for  both  geographic  and  celestial  target 
calculations  contain  the  geodetic  and  corrected  geomagnetic 
coordinates  of  the  vehicle's  position,  and  the  LVLH  azimuth 
and  elevation.  The  display  consists  of  the  azimuth  and 
elevation  versus  time,  with  the  angular  quantities  of  azimuth 
and  elevation  drawn  in  white  for  the  periods  that  the  earth 
did  not  obstruct  the  viewing  of  the  target  location.  For 
periods  of  obstruction  the  angles  are  drawn  in  red.  Each 
display  will  contain  the  length  of  time  selected  between  a 
maximum  of  12  hours  and  a  minimum  of  1  hour.  Periods  greater 
than  12  hours  are  shown  as  multiple  displays. 


6 . 4  AITS  Databases 


Various  databases  are  utilized  by  the  AITS  system  to  generate 
graphical  displays  and  determine  LOS,  FOV,  and  vehicle 
position.  The  users  create  their  own  set  of  parameters  by 
defining  attitude  angles  to  use,  auroral  contours,  gimbal 
orientation  angles,  and  time  periods  for  use  in  each  of  the 
analyses.  This  database  is  the  only  one  that  the  user  is 
allowed  to  alter  during  the  use  of  the  AITS  system  processes. 
All  other  databases  used  by  the  AITS  system  are  resident  on 
the  computer  system  as  a  master  directory  to  prevent 
inadvertent  alteration.  Sensor  specific  information 
concerning  the  sensor's  mount  angles  and  FOV  dimensions  are 
not  alterable  from  process  to  process  since  this  information 
is  specific  to  the  sensor. 


6.4.1  SAP  Stellar  Catalog 

The  Smithsonian  Astronomical  Observatory  ( SAO )  stellar 
catalog  consists  of  258,997  star  locations  referenced  for  the 
equinox  of  1950(7) .  Associated  with  each  data  location  is  a 
visual  magnitude  and  coefficients  to  correct  for  the  yearly 
motion  from  1950.  The  SAO  catalog  has  been  modified  to 
create  a  working  database  for  access  by  the  AITS  software 
package.  This  working  database  consists  of  almost  10,000 
stellar  locations  converted  to  the  year  and  month  of  launch. 

The  major  criterion  for  the  deletion  of  stellar  locations  in 
the  creation  of  the  working  database  was  the  inability  to  j 

view  these  locations  without  utilizing  specially  designed  j 

optical  sensors.  The  database  structure  created  requires 

minimal  computer  time  to  display  the  stellar  locations  that  I 

< 

are  within  the  FOV.  This  is  accomplished  by  the  use  of  a  ; 

direct  access  database  where  each  element  of  the  database 

defines  a  10  by  20  degree  grid.  1 

; 


6.4.2  World  Ma 


The  earth's  land  mass  representation  is  obtained  with  the  use 
of  two  separate  databases,  continental  outlines  and 
rasterized  panels  for  each  land  mass.  These  were  both 
created  from  a  database  maintained  by  the  Data  Systems  Branch 
at  the  Air  Force  Geophysics  Laboratory.  The  rasterized  world 
map  database  was  generated  from  the  original  by  separating 
the  database  into  grids  of  latitude  and  longitude  to  allow 
for  individual  grid  panel  rasterization.  This  database  is 
designed  as  a  direct  access  database  with  each  grid  panel 
being  accessed  as  a  separate  element  allowing  greater  access 
capability  by  the  camera  simulation  portion  of  the  AITS 
software . 


6.4.3  Planned  Mission  and  Updated  Kepler  Elements 

The  AITS  system  requires  availability  of  the  planned  mission 
state-vector  elements.  These  are  obtained  from  either  the 
Johnson  Space  Center  (JSC)  Super-Tape  microfiche  output  or 
from  the  Satellite  Control  Facility  ( SCF ) .  The  updating  of 
the  element  database  is  performed  on-orbit  by  the  user 
interactively  entering  the  actual  vehicle  state-vectors  at 
specific  times. 


6.4.4  Hardv/Gussenhoven  Auroral  Probabilitie 


This  database  consists  of  information  on  the  probability  of 
observing  high  flux  auroral  energy  for  a  single  orbital  pass 
of  the  auroral  region.  The  database  contains  two  categories 
of  probability  for  magnetic/solar  activity  level,  Kp  per  pas 
and  Equatorward  Boundary  ( EQ )  per  pass'1,2,31.  The  database 
was  obtained  from  the  analysis  of  sensor  data  derived  from 
DMSP  satellites.  The  difference  from  the  EQ  probability 
binning  was  the  use  of  midnight  auroral  boundary  in  place  of 
Kp. 

The  database  structure  consists  of  16  different  contours  for 
the  Kp  probability  and  25  different  contours  for  the  EQ 
probability.  Each  contour  represents  a  particular 
probability  condition  of  auroral  energy,  binned  by  the 
appropriate  criteria  at  either  two  or  three  percentage 
levels.  The  storage  format  allows  for  alteration  of  the 
database  given  new  contours  from  future  probability  study 
information. 


6.4.4. 1  Feldstein  Auroral  Oval 

This  database  consists  of  eight  auroral  oval  representations 
of  inner  and  outer  boundary  per  Feldstein  derived  auroral 
region'11.  Each  oval  is  stored  in  geomagnetic  coordinates 
of  local  time  and  latitude.  These  coordinates  are  converted 
and  displayed  by  the  Flight  Simulation  process. 


DMSP  Satellite  Images 


This  database  is  created  by  the  user,  during  on-orbit 
support,  from  Global  Weather  Center  (GWC)  communications 
consisting  of  DMSP  satellite  auroral  images.  Analysts  at  GWC 
evaluate  the  DMSP  auroral  image  and  transmit,  via  telephone, 
information  on  the  auroral  activity  region  viewed  by  DMSP 
satellites.  This  information  is  interactively  entered  into 
the  system  by  the  user  to  create  a  contour  database  for 
display  by  the  Flight  Simulation  process.  These  contours  are 
archived  on  the  AITS  computer  system  database  directory. 


6.4.6  Mission  Timeline 

This  information  is  entered  prior  to  launch  and  updated 
during  the  mission  as  required.  The  database  consists  of  all 
the  information  relevant  to  the  times  for  each  of  the  sensor 
measurements . 


6.4.7  Station  Coordinates 

Latitude  and  longitude  of  SCF,  GSTDN,  and  ADCOM  Radar 
Tracking  Stations  ( RTS  )  are  maintained  at  the  AITS  system 
master  directory.  These  locations  can  be  displayed  on  world 
maps  projections  if  required,  or  accessed  by  the  user  as 
ground  target  locations. 


6.5  Coordinate  Systems 


The  AITS  software  system  utilizes  a  number  of  coordinate 
system  in  the  generation  of  sensor  FOV,  LOS,  and  vehicle 
position.  The  most  used  systems  and  transformation  matrixes 
are  described  below. 


6.5.1  Local  Vertical  Local  Horizontal  ( LVLH ) 

This  coordinate  system  is  referenced  to  the  vehicle's 
position  and  direction  of  travel  (velocity  vector).  The 
primary  application  is  in  the  conversion  of  vehicle 
orientation  to  ECI  coodinates.  The  LVLH  coordinate  system  is 
a  right-handed  system  with  the  geocentric  radius  to  the 
vehicle  and  the  velocity  vector  defining  the  orbital  plane  at 
that  instant  in  time  (Figure  6.5).  The  LVLH  z-axis  lies 
along  the  geocentric  radius  vector  to  the  vehicle  and  is 
positive  toward  the  center  of  the  earth.  The  LVLH  y-axis  is 
the  normal  to  the  orbital  plane  and  is  determined  by  the 
cross-product  of  the  LVLH  z-axis  and  the  vehicle's  velocity 
vector.  The  cross-product  of  the  LVLH  y-axis  and  the  LVLH 
z-axis  generates  the  LVLH  x-axis,  completing  the  right-handed 
coordinate  system.  The  vehicle's  body  orientation  is 
transformed  by  the  LVLH  matrix  to  determine  the  vehicle's 
coordinates  in  an  ECI  coordinate  system. 


ll 


The  Euler  rotation  for  the  shuttle  attitude  is  a  pitch,  yaw, 
and  roll  sequence.  The  shuttle's  body  coordinate  system  is  a 
right  handed  system  with  the  x-axis  parallel  to  the  vehicle's 
structural  body  axis  (positive  toward  the  shuttle  nose), 
z-axis  parallel  to  the  vehicle's  plane  of  symmetry  and 
perpendicular  to  the  x-axis  (positive  down  with  respect  to 
the  shuttle  fuselage),  and  y-axis  completing  the  right-handed 
system  (positive  toward  the  shuttle's  right  wing).  The  Euler 
angles  for  the  shuttle  are  measured  positive  as  follows: 
pitch  up  (positive  x-axis  to  negative  z-axis),  yaw  right 
(positive  x-axis  to  positive  y-axis),  and  roll  clockwise 
(positive  y-axis  to  positive  z-axis)  (Figure  6.6). 


6.6  User's  Guide 


The  AITS  software  system  is  an  interactive  menu  driven 
package  that  allows  the  user  to  select  options  and  enter 
numeric  values  with  great  ease.  The  AITS  menus  require  a 
terminal  device  that  is  Digital  Video  Terminal  ( VT ) 
compatible.  The  AITS  system  utilizes  terminal  function  and 
cursor  movement  keys  in  the  selection  of  values  for  the 
processes.  The  menus  allow  the  user  to  either  use  the 
default  values  that  are  displayed  or  to  alter  the  value 
before  continuing  to  the  next  menu  option.  The  use  of  cursor 
movement  keys  allows  the  user  to  go  forward  or  backward  in 
the  menu  sequence  to  alter  or  verify  the  parameters  entered. 

The  AITS  system  is  initiated  by  the  user  entering  the  command 
" AFGL "  at  the  VAX  operating  system  prompt  "$"  after  login  on 
the  VAX.  The  main  menu  is  then  displayed  for  the  user  to 
select  the  desired  process  (Figure  6.7).  The  menu  display 
will  contain  the  AFGL  logo  if  the  terminal  is  a  Tektronix 
model  4109.  The  following  section  defines  the  processes  that 
the  user  may  select  from  the  main  menu. 
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1  -  UEHICLE  FLIGHT  SIMULATION 

2  -  EPHEMERIS  CALCULATIONS 

3  -  CELESTIAL  OPERATIONS  SIMULATION 

4  -  CAMERA  SIMULATION 

5  -  GIMBAL/OFF -TRACK  ANGLES 

6  -  SPECIAL  APPLICATIONS 
3  -  RETURN  TO  UAX'UMS 

PROCESS  DESIRED  - 


Figure  6.7  AITS  Main  Menu 


172 


References 


1.  James  A.  Mialen:  Auroral  Oval  Plotter  and  Ncmoyra^n  Determining  Corrected 
Geomagnetic  Local  Time,  Latitude,  and  Longitude  for  High  Latitudes  in  the 
Northern  Hemisphere.  Environmental  Research  Papers,  No.  327,  AFCRL-70-0422, 
July  1970,  AD713170. 

2.  W.J.  McNeil,  D. A.  Hardy,  R.R.  O'Neil:  Private  Communication. 

3.  W.J.  McNeil,  D.A.  Hardy,  R.R.  O'Neil:  Private  Communication. 

4.  M.S.  Gussenhoven,  D.A.  Hardy,  N.  Heineman:  Systematics  of  the  Equatorwara 

Diffuse  Auroral  boundary,  Journal  tor  Geophysical  Research,  88,  5692,  1983. 

5.  K.  Minka:  Orbit  Determination  and  Analysis  by  tne  Minimum  Variance  Method. 
Martin  Company,  Baltimore  Division,  Er  13950.  Prepared  tor  AFCRL,  OAR  (CRMXA) , 
USAE.  AECRL-65-579,  August  1965,  AD625453. 

6.  K.  Minka:  Orbit  Determination  and  Ephemeris  Computation.  Martin  Company, 
Baltimore  Division.  Prepared  for  AFCRL,  OAR  (CRMXA)  USAF.  AFCRL-66-259, 

May  1966,  AD637206. 

7.  Star  Cataloy:  Positions  and  Proper  Motion  of  258,997  Stars  tor  the  Epoch  and 
Equinox  of  1950.0.  Smithsonian  Institution,  Washington,  D.C. ,  1966. 


\ 


173 


Program  TTLOK,  an  adaptation  of  the  standard  AFGL  satellite 
ephemeris  program  LOKANGL,  provides  the  ability  to  compute 
ephemerides  for  the  double  thrust  situation.  A  double  thrust 
is  defined  as  two  successive  impulsive  thrusts  between  which 
no  intermediate  Keplerian  elements  or  P/V  (position/velocity) 
vectors  are  available.  Under  these  circumstances,  the 
standard  LOKANGL  program  is  unable  to  provide  valid  ephemeris 
calculations  for  the  time  interval  between  successive 
thrusts.  Experimenters'  need  for  coverage  during  this  period 
has  prompted  the  development  of  TTLOK  to  fill  the  void. 


Functional  Description 


The  reader  should  refer  to  Reference  1  for  an  overview  of  the 
basic  LOKANGL  program.  The  focus  of  activity  in  the  TTLOK 
modifications  to  LOKANGL  is  preparation  of  a  suitable  TAPE4 
file.  As  discussed  in  the  reference,  TAPE4  data  represents  a 
transformation  from  inputted  standard  Keplerian  element  sets 
or  P/V  (position/velocity)  vector  sets  to  a  related  set  of 
mean  Keplerian  elements  and  their  associated  time 
derivatives,  by  use  of  which  element  values  can  be  derived 
for  any  arbitrary  time  of  interest.  The  purpose  of  the  TTLOK 
modification  is  to  provide  the  information  needed  for 
spanning  the  time  interval  between  two  successive  thrusts 
during  which  no  orbital  measurements  are  available. 


Once  all  this  type  of  information  has  been  extracted  from  the 
input  (element  cards  or  P/V  vector  cards)  and  captured  in 
suitable  format  on  TAPE4,  TTLOK  (or  LOKANGL)  rewinds  the 


TAPE4  file  and  reads  the  first  record.  Each  record  consists 
of  a  set  of  elements,  their  derivatives,  their  epoch,  and  the 
latest  value  of  time  thru  which  these  data  are  to  be  applied 
to  ephemeris  calculations.  TAPE4  records  are  ordered 
chronologically  in  increasing  value  of  epoch  for  the 
associated  mean  element  set. 

For  each  instant  of  time  for  which  ephemeris  data  is 
requested,  beginning  with  the  earliest  and  proceeding 
sequentially  with  increasing  values  of  time,  LOKANGL  employs 
these  elements  and  their  temporal  derivatives  in  a  Taylor's 
series  expansion  to  calculate  mean  elements  for  each  of  the 
requested  instants  of  time.  The  instantaneous  values  of  mean 
elements  are  then  transformed  to  coordinate  values  defining 
the  vehicle's  instantaneous  location  and  velocity  (i.e.,  the 
ephemeris).  When  time,  the  independent  variable,  has 
advanced  to  the  end  of  the  region  of  applicability  of  the 
first  record  on  TAPE4 ,  the  next  record,  which  has  a  later 
cut-off  time  than  its  predecessor,  is  read  in.  The  forgoing 
process  then  continues  cyclically  until  the  entire  span  of 
time  for  which  calculations  have  been  requested  has  been 
covered.  (The  last  record  on  TAPE4  is  assigned  a  cut-off 
time  of  essentially  infinity,  permitting  calculations  to 
extend  as  far  as  desired  beyond  the  epoch  of  the  last  entry 
on  TAPE4. ) 

Standard  LOKANGL  is  unable  to  provide  ephemeris  values 
between  successive  thrusts  having  no  intermediate  point  of 
orbit  measurement.  TTLOK  eliminates  this  limitation  by 
providing  entries  on  TAPE4  which  provide  computational 
coverage  of  the  inter- thrust  period.  The  procedure  for 
generating  these  new  TAPE4  entries  involves  modeling  the 
thrust  process  as  follows: 


•  The  duration  of  the  thrust  is  assumed  infinitesimal. 
Thus  the  thrust  is  represented  as  an  acceleration 
impulse  occurring  at  the  nominal  thrust  time. 

•  The  thrust  impulse  (or  equivalently,  the  resultant 
discontinuous  change  in  vehicle  velocity)  is  assumed 
to  have  the  direction  of  the  vehicle  velocity 
immediately  before  thrusting. 


It  is  assumed  that  numerical  values  for  the  time  of  thrust 
and  the  associated  increment  in  velocity  are  available  from 
the  satellite  tracking  agency.  The  entries  on  TAPE4 
associated  with  the  double  thrusting  should  satisfy  the 
requirements  that: 

•  At  each  of  the  two  instants  at  which  thrusting  occurs 
there  should  be  discontinuous  changes  in  mean  elements. 

•  The  discontinuities  in  the  new  elements  should  reflect 
the  requisite  jump  in  value  of  the  vehicle  velocity  at 
the  time  of  thrust. 


•  Vehicle  position  should  remain  continuous  across  each  of 
the  thrust  times. 


7 . 3  Punched  Card  Changes  from  Standard  LOKANGL 

Input  cards  are  identical  to  those  for  standard  LOKANGL 
except  for  the  thrust  card,  to  which  is  added  the  increment 
in  velocity  in  columns  34-43  shown  below.  In  addition,  there 
are  now  two  successive  thrust  cards. 


Card 

No. 

Variable 

Name 

Card 

Col. 

Format 

Variable  Description 

2+  IYTH 

year  of  thrust 

19-20 

12 

Two  digit  value  of 

IDTH 

(or  IDTH2 ) 

21-23 

13 

Three  digit  value  of 

day  number  of  thrust 

SECTH 

(or  SECTH2 ) 

24-33 

F10.3 

Thrust  time  in 

seconds  of  day 

DELV1 

(orDELV2 ) 

34-43 

F10.3 

Increment  in  velocity 

in  feet  per  second 

Note  the  restriction  that,  as  presently  implemented,  TTLOK 
requires  that  the  sequence  of  input  element  (or  P/V)  set 
cards  includes  at  least  two  element  sets  prior  to  the  first 
thrust  card  and  at  least  two  element  sets  subsequent  to  the 
second  thrust  card.  This  requirement  is  not  a  fundamental 
limitation,  but  rather  is  a  consequence  of  the  particular 
design  approach  which  has  been  implemented. 


Ephemeris  output  from  TTLOK  is  identical  to  the  standard 
LOKANGL  print  options  (Reference  1).  The  header  has  been 
modified  to  accommodate  two  thrusts  and  to  indicate  the 
velocity  increment  for  each.  Elements  used  in  the 
inter-thrust  period,  and  associated  derivatives,  are 
printed  out  in  standard  LOKANGL  format. 
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7 . 5  Mathematical  or  Logical  Procedures 

As  noted  previously,  each  record  on  the  TAPE4  file  provides 
mean  Keplerian  elements,  their  derivatives,  the  epoch,  and 
DELDAX,  a  time  interval  which  is  to  be  added  to  the  epoch  to 
obtain  the  latest  time  for  which  data  in  the  record  is  to  be 
used  to  compute  an  ephemeris. 

a .  Background 

To  familiarize  the  reader  with  operation  of  the  LOKANGL 
system  several  simple  situations  will  be  illustrated. 

Consider  the  sequence  diagram  of  Figure  7.1.  There  is  only 
one  set  of  element  cards  in  the  input  deck.  They  are  denoted 
by  X  and  have  epoch  tx .  Correspondingly,  there  is  one  entry 
on  TAPE4  and  its  DELDAX  is  essentially  infinite.  Derivatives 
of  the  elements  are  estimated  based  on  geopotential  and  drag 
modules  internal  to  the  program.  This  single  TAPE4  entry  is 
used  for  calculations  for  all  times. 

Next  consider  the  case  of  two  element  sets  as  illustrated  in 
Figure  7.2.  DELDAX  for  the  first  entry  on  TAPE4  equals 
t2-t1.  The  derivatives  for  the  first  entry  represent  an 
interpolation  between  element  values  at  tx  and  those  at  t2. 

The  second  entry  on  TAPE4  has  an  infinite  DELDAX.  Thus  its 
usage  extends  from  t2  to  infinity.  The  derivatives  for  the 
second  TAPE4  record  are  obtained,  not  from  interpolation,  but 
from  the  internal  models. 

Consider  now  the  case  of  a  single  thrust  shown  in  the 
sequence  diagram  of  Figure  7.3.  Standard  element  cards  again 
represented  by  X,  and  the  thrust  card  by  T.  There  are  two 
element  sets  preceding  and  two  following  the  thrust.  The 
first  entry  on  TAPE4  has  tx/t2  interpolated  derivatives;  it 
is  used  between  minus  infinity  and  t=t2.  The  second  entry 
has  model-computed  derivatives  and  is  used  between  t2  and  t3. 
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gure  7.2  TAPE4  Structure  for  the  Case  of  Two  Element  Sets 


The  third  entry  has  derivatives  based  on  t4/t5  interpolation 
and  is  used  for  computations  between  t3  and  t5.  The  fourth 
entry  is  used  from  t5  onward  and  employs  model  derivatives. 

The  forgoing  examples  illustrate  the  general  rule  regarding 
applicability  of  element  sets  and  derivatives  contained  in  a 
given  TAPE4  record:  The  range  of  applicability  begins  at  the 
cut-off  point  of  the  preceding  record  (minus  infinity  for  the 
first  record  in  the  file)  and  extends  up  to  te+DELDAX  where 
te  is  the  epoch  for  the  given  record. 

b.  Double  Thrust  Outside  the  Inter-thrust  Interval 

The  sequence  diagram  for  the  case  of  a  double  thrust  is  shown 
in  Figure  7.4.  Standard  cards  with  epochs  t3  and  t2  are 
handled  in  the  usual  fashion  on  TAPE4:  i.e.,  a  t3  epoch  with 
t1/t2  derivatives  and  DELDAX  =  t2- tx .  But  when,  next,  TTLOK 
recognizes  that  the  t2  element  cards  are  followed  by  two 
successive  thrust  cards,  the  same  elements  and  derivatives 
are  written  a  second  time  to  TAPE4,  but  in  this  case  with 
DELDAX  =  t3- 13 .  This  repeated  record  written  to  TAPE4  provides 
coverage  within  the  pre-thrust  interval  from  t2  to  t3  using 
tj/t2  interpolated  derivatives.  Note  that  this  differs  from 
the  corresponding  situation  in  the  single  thrust  case  where 
model  derivatives  are  used  to  provide  coverage  between  the 
epoch  of  thrust  and  the  epoch  of  the  preceding  element  set. 

The  sequence  of  input  element  cards  depicted  in  the  figure 
results  in  5  records  being  written  to  TAPE4.  The  first  two, 
as  noted,  differ  only  in  their  values  of  DELDAX.  The  third 
record  covers  the  inter-thrust  period  and  is  discussed  below. 

Finally,  the  post-thrust  TAPE4  records  of  epochs  t5  and  t6 
are  similar  in  content  and  usage  to  the  two-card  case  of 
Figure  7.2.  The  only  difference  is  that,  in  the  present 
case,  usage  of  the  t5,t6  pair  begins,  not  at  minus  infinity, 
but  immediately  following  t4,  the  time  of  the  second  thrust. 


Structure  for  Double  Thrust 
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c .  Coverage  of  the  Inter-Thrust  interval 

The  procedures  for  obtaining  element  coverage  up  to  t3  (first 
thrust)  and  beyond  t4  (second  thrust)  have  been  described. 
Needed  now  is  coverage  within  the  interval  from  t3  to  t4 . 
t1/t2  extrapolated  elements  are  evaluated  for  the  time  just 
prior  to  t3  and  are  converted  to  osculating  P/V  vectors;  and 
similarly  for  the  time  just  following  t4 .  LOKANGL  routine 
MNTOPV  is  used  for  this  purpose.  Next,  it  is  assumed  that 
the  thrusts  create  increments  in  vehicle  velocity  which  are 
aligned  vectorially  with  the  pre-thrust  velocity  vector.  The 
V's  in  the  P/V  sets  at  t3  and  t4  are  incremented  (decremented) 
accordingly.  These  new  osculating  P/V  sets  are  then  trans¬ 
formed  by  LOKANGL  routine  PVTOMN  to  equivalent  mean  elements. 

Next,  derivatives  valid  in  the  t3  to  t4  interval  are 
calculated  using  the  standard  LOKANGL  code  for  evaluation  of 
derivatives  of  elements.  An  entry  is  then  written  to  TAPE4 
with  the  elements  being  the  post-thrust  #1  values  with  epoch 
t3;  derivatives  are  as  just  described;  and  DELDAX  =  t4-t3 .  The 
flow  of  the  overall  procedure  is  illustrated  in  Figure  7.5. 

The  procedure  for  handling  the  inter-thrust  period  is 
mechanized  in  subroutine  XFORM.  Given  a  set  of  mean 
elements,  their  epoch  (t),  their  derivatives,  a  prescribed 
time  interval  (At),  and  a  velocity  change  ( AV ) ,  XFORM 
performs  a  transformation  in  which  the  given  Kepler i an 
elements  are  converted  to  a  second  set,  at  epoch  t+At ,  t  ■ 
which  the  magnitude  of  the  vehicle  velocity  at  t  *  At  ha*- 
changed  by  an  amount  AV  from  the  value  which  would 
calculated  using  the  given  elements.  LOKANGL  nut  r 
MNTOPV  and  PVTOMN  are  employed  in  this  velocity 
transformation  process.  They  provide  the  me  • 
reversibly  converting  between  mean  Kepleii.c  . 
equivalent  osculating  P/V  vectors. 
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8.0  Software  for  Processing  Celestial  Aspect  Sensor  Data 

*  ■  . . ~  ’  "  _  ~  ~~  ‘  '  i 

( 

8 . 1  Introduction 

A  number  of  rocket-borne  experiments  require  high  accuracy 
attitude  information  to  successfully  interpret  in-flight 
measurements.  Instrument  line-of -sight  is  normally 
determined  from  gyroscopic  systems  mounted  in  the  vehicle 
payload.  However,  the  results  from  these  systems  may  contain 
errors  that  limit  the  quality  of  the  attitude  information. 

For  certain  demanding  applications,  gyro  system  measurements, 
alone,  fail  to  satisfy  the  requirements  for  precision.  In 
such  cases,  sun,  horizon,  and  celestial  aspect  sensors  can  be 
used  to  provide  correction  factors  to  be  applied  to 
gyroscopic  measurements  to  create  more  accurate  attitude 
profiles.  This  report  contains  a  general  description  of  the 
techniques  designed  and  implemented  to  calculate  corrected 
line  of  sight  information  for  a  payload  which  carried  a 
Celestial  Aspect  Sensor  (CAS). 

Utilization  of  the  CAS  measurements  required  the  development 
of  a  software  processing  system,  the  CAS  Processor  (CASP),  to 
convert  the  raw  CAS  data  into  corresponding  corrections  to 
the  nominal  vehicle  attitude  parameters.  Suitably  combined 
with  the  gyro  measurements,  these  CAS  updates  provide  an 
attitude  database  of  improved  accuracy.  This  upgraded 
attitude  profile  is  then  used  for  precise  calculations  of  the 
tangent  height  (and  related  geometrical  parameters) 
associated  with  the  field  of  view  of  the  prime  sensor,  an 
instrument  which  is  oriented  along  the  payload  body  axis. 


The  software  for  processing  of  CAS  data  has  been  developed  in 
two  stages.  The  first  of  these  is  documented  in  Reference  1, 
a  technical  memorandum.  The  second  stage,  which  is  the 
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subject  of  the  present  report,  represents  an  extension  of 
that  earlier  work  and  is  based  upon  a  data  processing 
algorithm  proposed  by  Dr.  H.A.  Miranda  and  Dr.  R.  Miranda, 
both  formerly  of  Epsilon  Co,  the  fabricators  of  the  CAS 
instrument.  The  effort  to  collect  and  process  CAS  data  was 
performed  in  support  of  the  Elias  rocket  experiment,  a  flight 
which  collected  IR  measurements  of  auroral  radiation. 

There  are  three  primary  types  of  data  upon  which  operation  of 
the  CASP  system  is  based.  In  addition  to  the  star  detections 
provided  by  the  CAS  instrument  itself,  the  other  databases 
are  the  vehicle  gyro  measurements  and  the  Smithsonian 
Astronomical  Observatory  (SAO)  star  catalog.  The  latter  is  a 
listing  (stored  on  magnetic  tape)  of  stars,  together  with 
their  magnitudes  and  locations  on  the  celestial  sphere  (i.e., 
right  ascension,  declination).  The  CASP  system  utilizes 
attitude  data  from  the  gyros  to  evaluate  the  approximate 
orientation  of  the  CAS  and  the  corresponding  field  of  view. 
The  star  catalog  data  can  then  be  used  to  identify  those 
stars  which  should  fall  within  the  field  of  view  of  the  CAS, 
together  with  their  locations  expressed  in  CAS  instrument 
coordinates.  Under  ideal  conditions  (no  instrument 
calibration  errors,  no  detection  false-alarms,  perfect  gyro 
data,  and  completeness  of  the  catalog)  all  of  the  instrument 
detections  should,  in  principle,  precisely  overlay  a  subset 
of  the  totality  of  catalog  stars  within  the  field  of  view. 

In  practice,  a  displacement  is  expected  between  the  spstial 
pattern  of  the  star  catalog  items  and  the  corresponding 
pattern  of  the  CAS  detections;  i.e.,  the  set  of  catalog 
items  and  the  set  of  detections  should  share  a  common  pattern 
of  geometrical  interrelationships  among  their  individual 
members,  but  there  is  a  displacement  between  these  two 
geometrically  similar  patterns.  For  purposes  of  data 
analysis,  it  is  assumed  that  such  displacements  are  due 
totally  to  gyro  measurement  error  (i.e.,  the  CAS  instrument, 
as  calibrated,  performs  perfectly).  The  objective  of  the 
analysis  performed  by  CASP  system  is  to  evaluate  attitude 


corrections  which  will  bring  the  catalog  data  into  "minimum 
offset"  agreement  with  the  detections.  Graphical  CRT 
displays  which  superimpose  both  the  SAO  catalog  items  and  the 
reported  detections  onto  the  CAS  field  of  view  are  used  to 
illustrate  the  effect  of  deviation  between  the  actual  vehicle 
attitude  and  the  attitude  deduced  from  (uncorrected)  gyro 
measurements.  An  example  is  shown  in  Figure  8.1. 

8.1.1  Description  of  Sensor 

The  Celestial  Aspect  Sensor  ( CAS )  is  a  scanning  instrument 
with  a  rectangular  field  of  view  of  approximately  21  by  25 
degrees.  The  instrument  uses  a  sawtooth  drive  to  scan  a 
mirror  up  and  down  the  21  degree  y-axis  of  the  field  of  view. 
As  the  mirror  scans,  it  reflects  the  stellar  emissions  into 
the  sensor.  The  optical  path  from  the  mirror  passes  through 
an  objective  lens  to  an  image  intensifier,  and  subsequently 
through  a  relay  lens  to  the  Charged  Coupled  Device  (CCD) 
detector.  The  CCD  detector,  a  256  element  linear  array  with 
an  S25  light  intensity  response,  creates  individual  raster 
lines  oriented  along  the  25  degree  x-axis  direction  of  the 
instrument.  The  mirror  requires  0.4  seconds  to  complete  a 
full  scan  and  thereby  to  create  a  256  by  256  pixel  field  of 
view  of  the  instrument.  Figure  8.2  is  a  cut-away  view  of  the 
construction  of  the  CAS. 

Every  0.1  seconds  the  on  board  processor  determines  the  x,y 
locations  of  the  eight  highest  intensity  pixel  responses 
located  in  the  one-fourth  of  a  field-of-view  swept  over  in 
this  partial  scan.  For  each  selected  detection  location,  the 
processor  then  calculates  three  quantities:  the  magnitude  of 
intensity,  delta-x,  and  delta-y.  Delta-x  and  delta-y  are 
position  corrections  for  detections  that  straddle  pixels,  and 
are  used  to  obtain  sub-pixel  location  accuracy  for  each 
detection. 


Figure  8.3  illustrates  a  typical  pixel-straddling  image.  The 


Figure  8.1  CRT  Display  of  Field  of  View  of  CAS  (bounded  by 
rectangular  frame)  Showing  Star  Catalog  Entries 
(+'s)  and  30  Detections  (at  centers  of  small 
rectangular  boxes). 
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following  relations  are  used  by  the  on-board  processor  to 
calculate  parameters  for  each  image: 

Delta-x  =  B„.  x 
Delta-y  =  BM_!  -  BMtl 

Magnitude  of  Intensity  =  (BM  +  B^)  -  (BMtl  +  BM,2) 

A  full  scan  field  of  view  consists  of  32  detection 
coordinates  created  from  four  consecutive  scan  segments  of 
the  same  direction  of  mirror  rotation.  Two  of  those  32  sets 
of  detection  coordinates  are  obtained  from  fixed  LEDs  (light 
emitting  diodes)  mounted  on  the  cover  lens  of  the  instrument. 
These  are  calibration  sources  used  to  adjust  the  pixel  field 
of  view  for  mirror  deviation  and  fluctuation  introduced  by 
the  mechanical  drive  system  of  the  mirror. 

Pulse  Code  Modulation  (PCM)  telemetry  is  used  to  transmit 
data  from  the  vehicle  to  the  ground  station.  Within  each 
telemetry  frame,  the  instrument  packs  and  transmits  data  for 
a  single  M-scan  segment  of  the  field-of-view.  The  frames 
occur  at  a  rate  of  one  every  0.1  seconds.  A  telemetry  frame 
consists  of  eight  sets  of  the  following  data:  detection  x  and 
y  coordinates,  magnitude  of  stellar  intensity,  delta-x,  and 
delta-y.  The  structure  of  the  telemetry  frame  is  illustrated 
in  Figure  8.4. 

8.1.2  Reduction  of  Haw  Data 

The  flight  data  was  recorded  on  instrumentation  tape  at  the 
rocket  range  and  sent  to  the  Telemetry  Data  Processing 
Section  of  AFGL  to  create  formatted  raw  sensor  data  tapes. 
These  are  unpacked,  and  their  quality  is  checked  to  create 
working  databases  for  use  by  data  reduction  software. 

The  tracking  data  is  processed  by  trajectory  determination 
software  to  produce  a  database  of  standard  position 
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parameters.  The  trajectory  data  and  the  raw  database  of 
gyroscopic  measurements  are  used  by  existing  attitude 
determination  software  to  produce  a  preliminary  flight 
history  of  attitude  defining  the  line  of  sight  for  each  of 
the  sensors  on  the  payload. 

The  SAO  star  catalog  is  unpacked,  and  each  entry  is 
calibrated  to  an  S25  light  intensity  response  range.  The 
epoch  for  the  flight  and  nominal  Celestial  Aspect  Sensor 
(CAS)  f ields-of-view  during  the  flight  are  used  to  define  a 
subset  of  the  full  SAO  catalog  known  as  the  SAO  working 
database.  The  CAS  working  database  is  calibrated  and 
adjusted  to  permit  reconstruction  of  fields  of  view  of  the 
CAS  during  the  flight. 

Figure  8.5  presents  an  overview  of  the  data  flow  supporting 
the  operation  of  the  CASP  software  system.  To  initiate 
processing  of  CAS  data,  the  CASP  software  system  extracts, 
from  the  CAS  working  database,  four  consecutive  telemetry 
frames  of  the  same  scan  direction  to  reconstruct  a  CAS  field 
of  view  of  a  segment  of  the  celestial  sphere. 

Note  that  delta-x  and  delta-y  corrections  have  already  been 
applied  to  the  detection  readouts  by  an  on  board  processor  to 
yield  interpolated  coordinates  for  detections  that  straddle 
adjacent  pixels;  equations  1,  2,  and  3  of  Table  8.1 
illustrate  how  these  adjustments  are  applied. 

Each  field  is  corrected  for  mirror  misalignment  using  the  two 
reference  star  coordinates  contained  in  the  field  of  view 
(Eq.  4,  5).  The  remaining  readouts  are  calibrated  into 
instrument  degrees  and  corrected  for  magnification  distortion 
in  the  x-axis  direction  and  mirror  mounting  angle  skewness  in 
the  y-axis  direction  (Eq.  6,  7,  8). 
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Table  8.1 

Calibration  Equations 


Eq.  1 


Eq.  2 


Eq .  3 


Eq.  4 


If  byte( 6 )  <=  128: 

CM  =  byte ( 6 ) 

Else: 

CM  =  byte( 6 )  -  256 

SM  =  byte(l)  +  ( 256 ) ( byte( 2 ) ) 

If  byte(6)  <=  128: 

Cdelta-y  =  [byte( 6 )/( 2( byte  (1,2))  + 

byte ( 6 ) / ( 4 ( byte (1,2)  -  byte(6)))] 

Else: 

Cdelta-y  =  [ (byte( 6 )-256 )/( 4( byte( 1, 2 )  - 

(CM)))  +  ( byte( 6 )-256 ) /( 2( byte( 1,2) 
-  2 ( CM )  ) )  ] 

If  direction  bit  =  1  :  Sign  =  1 
Else:  Sign  =  -1 

Correct-Y'  =  byte(4)  +  Sign( Cdelta-y ) 

Correct-X'  =  byte(3)  -  f ( byte( 5 )/SM ) 
where  f  is  the  function  defined  in  Table  8.2 

Refl-X  =  X-coordinate  of  LED  #1 
Ref2-X  =  X-coordinate  of  LED  #2 
Refl-Y  =  Y-coordinate  of  LED  #1 
Ref2-Y  =  Y-coordinate  of  LED  #2 
Adj-X  =  (Refl-X  +  Ref 2-X ) /2  -  208 
Adj-Y  =  (Refl-Y  +  Ref2-Y)/2  -  113 

Correct-X''  =  Correct-X'  -  Adj-X 
Correct-Y' '  =  Correct-Y'  -  Adj-Y 
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Table  8.1 
( continued ) 

Eq.  6.  Correct-X  =  Correct-X'  -  1 . 863E-6( | Correct-X '  - 
140 | 3 ) 

Eq.  7.  Correct-Y  =  Correct-Y' '  -  0.01375 (Correct-X'  - 
140) 

Eq.  8.  Deg-X  =  Correct-X( 0. 001719 )180/pi 

Deg-Y  =  (255  -  Correct-Y )(0. 001417 ) 180/pi 
Eq.  9.  LOSDEC  =  declination  from  gyro  calculation 

LOSRA  =  right  ascension  from  gyro  calculation 
Fi  =  tan"1  (Deg-X  -  13.2992772902)/ 

(Deg-Y  -  13.9406770361)) 

Rads  =  (Deg-X  -  13 . 2992772902 )cos( Fi ) 

Eq.  10.  Inst-DEC  «  LOSDEC  +  (Rads)sin(Fi  -  Alpha) 

Eq.  11.  Inst-RA  =  LOSRA  +  (Rads)cos(Fi  -  Alpha) 


Table  8.2 

Discrete  Points  Defining  Function  f 
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8.2  Processing  Software  For  The  CAS 


A  celestial  aspect  sensor  (CAS)  was  flown  on  the  Elias  rocket 
flight  to  provide  high  precision  attitude  data  to  augment  the 
data  collected  using  horizon  sensor  and  gyro  instruments.  At 
the  most  rudimentary  level,  the  CAS  may  be  described  as 
consisting  of  a  linear  array  of  256  CCD  elements  together 
with  a  rotating  (actually  rocking)  mirror  which  provides  a 
scanning  capability.  The  period  of  the  scan  is  0.4  seconds, 
and  the  scan  direction  is  reversed  on  successive  frames. 

Each  scan  provides  a  two-dimensional  image  of  a  section  of 
the  celestial  sphere.  It  is  assumed  that  the  CAS  and  the 
rocket  body  jointly  constitute  a  single,  perfectly  rigid  body 
of  precisely  known  geometry.  Thus,  if  the  orientation  of  the 
CAS  is  known,  the  attitude  of  the  vehicle  is  deducible;  and 
conversely.  Further,  it  is  assumed  that  at  each  instant  of 
time  of  interest,  the  attitude  data  deduced  from  gyro 
measurements  is  sufficiently  accurate  to  ensure  that 
calculated  expected  locations  of  stars  within  the  viewing 
field  of  the  CAS  are  of  "ballpark"  caliber.  Basically,  the 
processing  system  evaluates  these  expected  star  locations, 
forms  pair-wise  associations  between  detections  and  catalog 
stars,  evaluates  the  displacements  between  the  observed  star 
detections  and  their  respective  star  catalog  counterparts, 
and  calculates  the  increments  in  attitude  angles  needed  to 
bring  the  two  sets  of  points  as  close  to  overall  coincidence 
as  possible.  Thus,  it  is  tacitly  assumed  that  mismatch 
between  calculated  star  locations  and  the  corresponding 
detection  points  is  caused  solely  by  error  in  the  gyro 
measurements.  Hence,  the  increments  in  attitude  angles 
needed  to  achieve  optimum  coincidence  can  be  viewed  as 
corrective  updates  to  the  gyro  data  base,  valid  at  the  time 
corresponding  to  that  of  the  CAS  scan.  (Actually  scanning  is 
not  an  instantaneous  process,  but  occurs  over  a  finite  period 
of  time,  during  which  the  payload  attitude  may  vary.  See 
Section  8.2.2  where  this  effect  is  taken  into  account.) 


The  operation  of  CASP  is  based  on  a  frame-by-frame  analysis 
of  the  CAS  data,  coupled  with  feed-forward  updating. 

Attitude  data  deduced  from  gyro  measurements,  alone,  is  never 
of  sufficient  accuracy  to  yield  perfect  alignment  between 
star  catalog  data  and  the  corresponding  CAS  detections. 
However,  if  the  correct  attitude  were  somehow  known  for  frame 
N,  then  that  attitude  data  could  be  updated  to  the  time  of 
frame  N+l  using  available  rate  gyro  data  in  a  Taylor's  series 
expansion.  It  is  assumed  that  attitude  data  obtained  by 
updating  from  an  adjacent  frame  is  of  sufficient  accuracy  to 
yield  near-coincidence  between  CAS  detections  and  their 
corresponding  catalog  items  (if  any).  Reported  detections 
failing  to  have  a  nearby  catalog  counterpart  can  be  assumed 
to  be  noise  and  can  be  ignored.  Thus  a  small  acceptance  cell 
can  be  erected  about  each  reported  detection.  (These  are  the 
30  rectangules  evident  in  Figure  8.1.  Each  detection  can 
then  be  matched  to  that  star  (if  any)  falling  within  its 
acceptance  cell.  If  more  than  one  catalog  item  falls  within 
a  given  cell,  the  most  probable  is  selected  (e.g.,  on  the 
basis  of  star  magnitude).  Small  attitude  adjustments  are 
then  evaluated  which  minimize  the  overall  misalignment  of 
detections  and  stars.  This  completes  the  analysis  of  frame 
N+l.  The  process  is  then  repeated  for  frame  N+2. 

The  forgoing  process  is  the  basis  for  the  automatic  mode  of 
operation  of  the  CASP  system.  Figure  8.6  illustrates  the 
operation  of  the  automatic  mode.  Once  initiated,  this  mode 
is  designed  to  process  successive  frames  sequentially, 
provided  only  that  the  quality  of  the  raw  CAS  data  is 
adequate  to  sustain  the  operation.  The  caliber  of  the  data 
can  be  judged  on  the  basis  of  the  percentage  of  the  30 
detections  reported  per  frame  which  represent  legitimate  star 
detections.  If  their  number  is  too  few,  the  system  will  have 
only  noise  to  which  to  lock  and  the  results  will  be  erratic 


Lgure  8.6:  Highly  Schematic  Representation  of  Normal  Initiation 
and  Subsequent  Automatic  Mode  of  CASP  Operation 


and  erroneous.  Internal  tests  have  been  provided  to  identify 
this  situation  and  to  terminate  automatic  processing  when  it 
occurs. 

As  noted  above,  gyro  data  alone  is  too  inaccurate  to  ensure 

near-coincidence  of  individual  detections  with  their 

corresponding  catalog  items.  How  then  can  the  automatic  mode 

be  initiated?  The  solution  is  a  manual  designation  process. 

The  analyst  is  provided  with  a  CRT  graphic  display  showing 

both  the  catalog  stars  and  the  detections  plotted  on  the  CAS 

field  of  view,  as  shown  in  Figure  8.1.  The  operator  judges 

which  detections  are  real  and,  using  a  cursor,  manually 

associates  them  with  the  corresponding  stars.  CASP  then 

evaluates  the  attitude  corrections  needed  to  minimize  the 
« 

overall  star/detection  offset,  and  automatic  processing  can 
then  proceed  unaided.  There  are  two  basic  differences 
between  the  manual  and  the  automatic  mode.  The  first  relates 
to  the  gyro  data  used  in  coordinate  transformations.  For  the 
manual  mode,  only  raw  gyro  data  is  employed.  In  the 
automatic  mode,  the  attitude  data  employed  is  an  amalgam  of 
the  raw  gyro  data  together  with  corrective  upgrades  obtained 
from  the  processing  of  a  time-adjacent  frame.  The  second 
difference  relates  to  the  process  of  associating  detections 
with  star  catalog  items.  In  the  manual  mode,  the  operator 
designates  matching  pairs  using  a  graphic  display  together 
with  a  cursor.  In  the  automatic  mode,  the  process  is 
automatically  performed  using  a  criterion  based  upon  the 
acceptance  cell  concept. 

The  following  subsections  provide  more  detailed  background  on 
the  operation  of  CASP. 

8.2.2  Software  System  Overview 


Figure  8.7  illustrates  the  major  functions  performed  in  the 
processing  of  a  single  frame  of  CAS  data.  The  process 
involves  transforming  the  coordinates  of  each  candidate  star 


from  the  right  ascension/declination  inertial  system  to 
instrument  coordinates  (x,y)  appropriate  for  the  given  time. 
This  transformation  is  a  function  of  the  attitude  of  the 
vehicle  at  the  given  time.  It  should  be  recognized  that  0.4 
second  of  time  elapses  during  the  collection  of  one  frame  of 
data.  This  is  the  time  required  for  the  mirror  to  complete 
the  y-direction  scan.  During  this  period  the  vehicle 
attitude  can  vary  appreciably.  In  order  to  identify  what 
attitude  values  apply  for  a  given  catalog  star,  its  time  of 
"detection"  must  be  evaluated.  But  time  within  a  scan  can  be 
measured  in  proportion  to  the  amount  of  mirror  rotation 
experienced:  i.e.,  in  terms  of  the  y  coordinate. 

Inter-fr.  tie  motion  is  treated  by  performing  a  two-pass 
coordinate  transformation.  On  the  first  pass,  attitude  data 
corresponding  to  the  start-time  of  the  frame  is  used.  This 
defines  an  approximate  y-coordinate .  That  value  of  y,  being 
obtained  through  the  process  of  mirror  rotation,  defines  the 
approximate  proportion  of  the  0.4  sec  scan  period  elapsed  at 
the  instant  when  the  star  would  be  detectable  by  the  CAS. 

With  this  new  time-reference  available,  the  attitude  values 
for  use  with  each  individual  star  are  updated  to  the  "time  of 
detection" .  The  coordinate  transformation  for  each  catalog 
star  is  then  performed  a  second  time  using  as  input  the 
attitude  appropriate  for  the  time  of  detection.  The  result 
is  a  set  of  x,y  coordinates  for  each  star  which  have  been 
corrected  for  inter-frame  motion.  In  principle,  the  process 
could  be  repeated  iteratively,  with  successively  more  exact 
results.  In  practice,  though,  one  correction  is  sufficient. 

The  process  of  calculating  an  attitude  correction  from  a 
given  frame  of  data  begins  with  centering  an  acceptance  box 
on  each  of  the  30  star  detections  reported  for  the  frame. 

Each  catalog  star  is  examined  to  determine  whether  it  falls 
within  any  of  these  boxes.  Box  occupancy  serves  to  define 
the  association  between  the  catalog  star  and  the 


corresponding  detection. 


Once  all  associations  have  been  identified,  a  mismatch  error 
function  is  defined  as  shown  in  Figure  8.7.  This  function 
serves  as  an  objective  function  for  the  subsequent 
optimization  (minimization)  process.  The  basic  form  of  the 
error  is  given  by  the  expression 


Mismatch  Error  = 


l  [WRi(X.tar  -  Xdat)?  + 


Wyi(Y,tar  -  Ydet)?] 


where  the  index  i  identifies  individual  associated 
star/detection  pairs,  the  summation  is  over  all  such  pairs, 
(x.t.r  -  xdet)i  and  ( yatar  ~  Ydet>i  are  the  x,y 
star/detection  coordinate  differences  for  the  i-th 
star /detection  association,  and  the  W's  are  weighting 
factors.  The  mismatch  error  is  functionally  dependent  on  the 
vehicle  attitude  parameters.  For  a  perfect  CAS  instrument 
and  a  set  of  observed  detections  that  represent  only  valid 
star  "hits"  and  no  noise-induced  "false  alarms",  the  mismatch 
error  function  reduces  (in  principle)  to  a  minimum  value  of 
zero  when  evaluated  for  those  attitude  parameters  which 
represent  the  true  orientation  of  the  vehicle. 


The  final  step  in  processing  a  CAS  frame  of  data  is  a  formal 
optimization  procedure  which  yields  incremental  corrections 
to  the  pitch,  yaw,  and  roll  attitude  parameters,  jointly 
evaluated  to  minimize  the  total  mismatch  error. 


8.2.3  Transformations  between  Coordinate  Systems 

Stars  are  organized  and  identified  in  the  Smithsonian  Catalog 
by  their  ECI-based  right  ascension  and  elevation.  However 
their  physical  manifestations  are  observed  by  the  CAS  in  the 
instrument's  f ield-of-view  coordinates.  Correlating  those 
two  bodies  of  data  entails  conversion  to  a  common  coordinate 


system.  The  major  instrument  of  the  ELIAS  payload  is  an  IR 
sensor  with  its  f ield-of-view  oriented  in  the  direction  of 
the  rocket  axis.  This  sensor  collects  data  on  the  IR 
emissions  of  the  auroral  regions  toward  which  it  is  directed. 
Consequently,  a  key  input  to  analysis  of  the  experimental 
data  is  identification,  as  a  function  of  flight  time,  of  the 
volume  of  space  under  observation.  This  viewing  region  at 
any  given  time  is  determined  by  the  instantaneous  orientation 
of  the  rocket  body  coordinate  system,  which  can  be  calculated 
from  attitude  measurements,  referenced  to  the  launcher. 
Ground-based  radar  observations  of  the  payload  provide  its 
location  expressed  in  coordinates  centered  at  the  radar  site. 

Payload  attitude  and  position  data  can  be  combined  to  yield 
the  tangent  height  data  which  is  needed  to  identify  the 
region  under  observation.  The  geometry  is  illustrated  in 
Figure  8.8.  Tangent  height  is  defined  as  the  minimum 
distance  between  the  rocket  line-of-sight  (LOS)  and  the 
surface  of  the  oblate  earth.  It  is  specified  in  terms  of 
that  minimum  distance  together  with  the  longitude  and 
geodetic  latitude  of  the  corresponding  point  on  the  surface. 
Here  again,  the  required  geometric  analysis  utilizes  multiple 
coordinate  system  conversions. 


As  the  forgoing  discussion  illustrates,  the  processing  of 
the  CAS  measurements  requires  that  the  software  system 
provide  the  capability  to  transform  freely  among  a  variety  of 
coordinate  systems.  The  following  subsections  describe  the 
major  transformations  implemented  in  the  CASP  software. 


8 . 2 . 3 . 1  Earth-Centered-Inertlal  to  Inertial -Fixed-Geographic 


The  earth-centered  inertial  (ECI)  system  is  defined  by  an 
x-axis  aligned  with  the  vernal  equinox  and  a  z-axis  oriented 
northward  along  the  earth's  axis  of  rotation.  Like  the  ECI, 
the  inertially  fixed  geographic  (IFG)  system  is  earth 
centered,  and  employs  the  earth's  rotational  axis  as  its 
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Figure  8-8  Illustration  of  Tangent  Height 


z-axis;  but  its  x-axis  is  positioned  in  the  Greenwich 
meridian  plane  at  the  instant  of  gyro  uncage  (taken  to  be 
launch  time). 


Note  that  because  it  is  defined  to  be  inertial,  the  IFG 
system  does  not  rotate  with  the  earth.  It  coincides  with  the 
conventional  rotating  geographic  coordinate  system  only  at 
gyro  uncage  time.  Figure  8.9  illustrates  the  relationship 
between  the  ECI  and  the  IFG  systems. 

8. 2. 3. 2  IFG  to  Vertical-East-North  ( VEN ) 

The  VEN  system  is  inertially  fixed,  is  centered  at  the 
launcher’s  location  ( topocentric ) ,  and  has  its  axes  pointing, 
respectively: 

°  vertically  (i.e.,  direction  of  local  plumb  line  on  the 
oblate  earth), 

°  east,  and 

°  north 

Like  the  IFG  system,  the  VEN  system,  being  inertially  fixed, 
does  not  participate  in  the  earth's  rotation;  it,  too,  is 
defined  in  terms  of  the  value  of  the  Greenwich  hour  angle  at 
the  time  of  uncage.  Figure  8.10  illustrates  the  relationship 
between  the  IFG  and  the  VEN  systems.  Note  that  account  must 
be  taken  of  the  oblateness  of  the  earth  in  specifying  the 
center  of  the  VEN  system. 

8. 2. 3. 3  VEN  to  Launcher  System 

The  VEN  and  Launcher  systems  share  a  common  center.  The 
Launcher  system  is  inertially  fixed,  with  axes  coincident 
with  the  rocket  body  axes  at  the  uncage  time.  Thus  the 
launcher  system  is  identical  with  the  body  system  at  launch 
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Figure  8.10  Geometry  Relating  the  IFG  and  the  VEN  Coordinate  Systems 
(V  is  perpendicular  to  surface  of  oblate  earth) 


time.  Subsequently,  the  rocket  body  system  undergoes 
rotational  and  translational  motion  associated  with  flight. 
The  launcher  system,  on  the  other  hand,  remains  inertial ly 
fixed  and  provides  the  reference  directions  against  which  the 
gyro  system  measures  in-flight  variation  in  rocket  body 
attitude. 

The  body  system  is  defined  as  follows: 

°  x-axis  is  the  roll-axis 

°  y-axis  is  the  pitch-axis 

°  z-axis  is  the  yaw-axis 

Figure  8.11  illustrates  the  geometry. 

8. 2.3.4  Launcher  to  Rocket  Body  System 


The  orientation  of  a  rigid  body  can  be  defined  in  various 
ways:  direction  cosines  (9  quantities);  quaternions  (4 
quantities);  Euler  angles  (3  quantities);  and  the  pitch,  yaw, 
and  roll  attitude  angles.  Common  rocketry  practice  is  to  use 
the  attitude  angles,  which  correspond  to  the  minimum  number 
of  independent  parameters  required  to  uniquely  specify 
attitude.  Each  of  the  pitch,  yaw,  and  roll,  angles 
represents  a  rotation  about  one  of  three  orthogonal  axes. 
These  are  finite  (i.e.,  not  infinitesimal)  angle  rotations; 
and,  thus,  they  do  not  commute. 

In  effect,  the  general  transformation  matrix,  which  relates 
an  arbitrary  instantaneous  attitude  of  a  rigid  body  to  a 
prescribed  reference  coordinate  system  (in  our  case,  the 
launcher),  can  be  decomposed  into  the  product  of  three 
matrices,  each  of  which  represents  a  rotation  about  an  axis 
through  prescribed  angles  (the  pitch,  yaw,  and  roll  angles). 
However,  tacitly  associated  with  a  given  decomposition. 


represented  as  numerical  values  of  pitch,  yaw,  and  roll,  is  a 
prescribed  sequence  in  which  the  corresponding  axial 
rotations  must  be  applied  in  order  to  produce  the  intended 
body  orientation.  The  non-commuting  of  these  three  rotations 
implies  that  the  correct  orientation  of  a  body  can  be 
obtained  from  given  values  of  pitch,  yaw,  and  roll  only  if 
the  corresponding  axial  rotations  are  applied  in  the  proper 
sequence.  For  the  Elias  flight,  that  sequence  is  pitch, 
first,  then  yaw,  and  finally  roll.  Figure  8.12  presents  the 
relationship  between  the  launcher  and  the  body  systems. 

8. 2. 3. 5  Rocket  Body  Coordinates  to  Instrument  x,y  Coordinates 

It  will  be  recalled  that  the  CAS  instrument  is  comprised  of  a 
linear  array  of  CCD  elements,  together  with  a  mirror  which 
provides  scanning  in  the  transverse  direction.  For  each  of 
the  30  detections  declared  per  frame,  x  and  y  coordinates  are 
provided.  Calibrated  in  terms  of  angle,  these  numbers 
identify  an  apparent  direction  of  look  associated  with  the 
corresponding  detection.  The  x,y  coordinates  are  referenced 
to  the  CAS  boresight  and  represent  an  angular  offset  from 
that  direction.  The  CAS  was  mounted  on  the  rocket  body  such 
that  its  boresight  direction  was  in  the  plane  of  the  body's 
roll  and  yaw  axes  with  the  x-direction  perpendicular  to  the 
roll-yaw  plane.  The  boresight  was  offset  by  an  angle  « 
(nominally  45°)  from  the  roll  axis.  The  mounting  arrangement 
is  illustrated  in  Figure  8.13. 


To  complete  the  chain  of  transformations  between  the  optical 
measurements  and  the  corresponding  star  locations  defined  in 
an  inertial  system,  we  need  the  transformation  between  x,y 
instrument  space  and  the  rocket  body  coordinate  system. 

Figure  8.14  illustrates  the  geometry.  A  point  p(x,y) 
represents  a  detection  defined  in  terms  of  its  instrument 
coordinates  x,y.  Contours  of  constant  x  are  circles  whose 
planes  parallel  to  the  plane  formed  by  the  yaw  and  roll  axes. 
Contours  of  constant  y  are  great  circles,  the  plane 


Figure  8-11  The  Body  Coordinate  System 


Figure 


-12  Relationship  between  the  Local  VEN  System  and  the  Launcher  Coordinate 
System  (X,Y,Z) 


of  each  of  which  contains  the  body  pitch  axis,  P.  The  point 
p  also  has  the  body  coordinates  Rp,Pp,  and  Yp.  The 
relationship  between  the  x,y  system  and  the  R,P,Y  system  is 
given  by 

sin  x  =  Pp 

tan  y  =  Yp/Rp 

8.2.4  Minimum  Displacement  Star-to-Detection  Fitting:  An 
Optimization  Process 

The  crucial  step  in  utilizing  CAS  observations  to  provide 
corrective  updates  to  gyro-based  attitude  data  is  to  evaluate 
those  incremental  attitude  changes  which  optimally  bring  the 
coordinates  of  the  CAS  detections  into  (near)  coincidence 
with  the  coordinates  of  the  associated  catalog  stars.  The 
independent  variables  in  the  process  are  the  three  attitude 
parameters:  pitch,  yaw,  and  roll.  The  analysis  is  performed 
in  the  f ield-of-view  coordinates  of  the  CAS  instrument. 

The  procedure  is  to  define  an  analytic  "objective"  function 
which  quantitatively  measures  the  degree  of  mismatch.  In 
CASP  the  objective  function  is  the  sum  of  the  squares  of  the 
individual  displacement  distances  in  the  x-y  plane  of  the  CAS 
coordinate  space.  The  Fletcher-Powell  optimization  algorithm 
(Reference  2)  is  then  employed  to  search  iteratively  for  that 
set  of  attitude  parameters  which  jointly  minimize  the 
objective  function.  This  algorithm  belongs  to  the  gradient 
search  class  of  optimization  procedures.  When  the  function 
to  be  optimized  is  ill-behaved,  problems  can  be  experienced 
with  this  type  of  approach.  However,  no  difficulty  has  been 
encountered  in  the  CASP  application. 

The  user  must  supply  the  Fletcher-Powell  algorithm  with  an 
expression  for  the  objective  function,  known  as  FUNCT,  by 
means  of  which  both  the  function  value  and  the  components  of 


the  gradient  can  be  evaluated.  Figure  8.15  illustrates  the 
use  of  FUNCT  to  compute  the  data  needed  to  support  rapid 
descent  minimization  algorithm. 

8.2.5  Tangent  Height  Calculation 

The  geometry  of  the  problem  to  be  addressed  is  illustrated  in 
Figure  8.8.  The  vehicle,  located  at  point  V  at  a  radial 
distance  of  OV  from  the  center  of  the  earth,  0,  has  its  roll 
axis  oriented  in  the  direction  R.  The  problem  is  to  locate 
the  point  of  closest  approach  between  the  surface  of  the 
oblate  earth  and  the  straight  line  through  point  V  in  the 
direction  of  the  vehicle's  center  line,  R.  At  the  point  of 
closest  approach  the  perpendicular  distance  between  this  line 
and  the  earth's  surface  is  defined  as  the  tangent  height.  In 
the  drawing,  the  Z  axis  is  the  earth's  rotational  axis  and 
the  x-y  plane  is  the  earth's  equatorial  plane. 

With  the  coordinates  of  point  V  and  the  components  of  R 
specified,  points  along  the  line  through  V  in  the  direction  R 
are  a  one  parameter  family,  that  parameter,  x,  being  distance 
along  the  line.  Similarly,  if  the  surface  of  the  earth  is 
modeled  as  an  oblate  spheroid,  z  coordinates  of  points  on 
that  surface  can  be  represented  as  functions  of  their  x-y 
coordinates.  An  expression  can  then  be  written  for  the 
distance  between  an  arbitrary  point  on  the  line  and  an 
arbitrary  point  on  the  surface.  This  function  depends  on 
three  independent  variables:  x,  y,  and  distance  along  the 
line.  The  method  of  solution  is  to  employ  this  function  as 
an  objective  function  which  is  to  be  minimized  using  the 
Fletcher-Powell  method  of  rapid  descent. 

The  minimization  process  is  initiated  from  the  spherical 
earth  solution  represented  by  the  radial  from  the  center  of 
the  earth  which  is  the  perpendicular  to  the  tangent  line 
(i.e.,  R  extended).  The  use  of  a  good  first  approximation  to 
initialize  the  process  tends  to  reduce  the  number  iterations 


mplementation  of  Subroutine  FUNCT.  With  suitably  selected  value 
>P»  AY,  and  aR  (incremental  values  of  attitude  parameters),  FUNCT 
s  used  to  provide  numerical  values  of  the  objective  function  and 
ts  derivatives. 


needed  to  achieve  a  solution. 


Quantities  computed  are  the  tangent  height,  Ht  and  geodetic 
and  geocentric  coordinates  of  the  point  on  the  spheroidal 
earth  below  the  point  of  tangency. 

8.2.6  Data  Quality  Considerations 

The  process  of  detecting  stellar  radiation  is  one  which  is 
inherently  less  than  perfect.  Viewed  simplistically,  the 
level  of  performance  is  a  function  of  factors  such  as  the 
intensity  of  radiation  of  the  individual  stars,  the 
background  "noise"  level  (i.e.,  the  dark-current)  of  the 
detectors,  the  effective  integration  time,  and  the  level  used 
as  a  detection  threshold.  Imperfection  in  the  detection 
process  can  manifest  itself  in  several  ways:  e.g., 

a)  false  alarms  (noise  only  events  incorrectly  declared 
to  be  detections) 

b)  missed  detections 

c )  incorrect  evaluation  of  the  coordinates  for  actual 
and/or  apparent  pixel-straddling  events 

The  CASP  software  system  performs  its  functions  too  far 
downstream  from  the  detection  process  to  be  able  to  affect 
either  b)  or  c).  However,  the  fact  that  the  telemetered  CAS 
data  involves  multiple  detection  events,  coupled  with  certain 
necessary  conditions  which  must  be  met  by  physically  real 
data,  afford  the  possibility  for  the  CASP  software  to  provide 
some  discrimination  against  false  alarms.  Features  in  the 
software  which  capitalize  on  these  principles  are  discussed 
in  the  following  four  subsections. 

8. 2. 6.1  Roll  Rate  Discrimination 

i 

i 


The  possibility  exists  for  the  system  to  lock  to  real  stars 
on  a  given  frame,  but  on  the  succeeding  frame  to  lock  to  one 
or  more  false  alarm  events,  yielding  erroneous  attitude 
updates  for  the  second  frame.  One  consequence  of  this 
situation  would  be  the  likelihood  that  the  magnitude  of  the 
change  in  roll  within  the  time  interval  between  the  first  and 
second  frames  would  exceed  the  bounds  imposed  by  the  laws  of 
dynamics  governing  motion  of  the  vehicle. 

A  test  is  performed  by  CASP  to  identify  excessive  inter-frame 
roll  increments.  When  such  an  occurrence  is  detected, 
automatic  frame-to- frame  processing  is  terminated;  and  manual 
initiation  is  needed  to  re-start  the  automatic  mode  of 
operation. 

8. 2. 6. 2  Use  of  Weighting  within  the  Objective  Function 

Consider  the  analysis  of  data  for  a  frame  for  which  n 
star-detection  associations  have  been  identified.  Among  the 
n  detections  there  is  a  range  of  detected  signal  levels. 
Likewise,  there  will  be  variation  in  optical  intensity  among 
the  n  stars.  A  crude  measure  of  probability  of  physical 
reality  of  associated  pairs  is  provided  by  a  table  such  as 
the  following: 


Probability  Event 
Star  Intensity  Detection  Level _ Is  Real _ 


High 

High 

High 

High 

Low 

Intermediate 

Low 

High 

Intermediate 

Low 

Low 

Low 

The  objective  function  which  is  to  be  minimized  in  order  to 
evaluate  the  "true"  vehicle  attitude  is  the  sum  of  the 
squares  of  the  individual  mismatch  "distances"  in  x-y  space 


for  each  of  the  identified  star-detection  associations.  To 
minimize  the  contributions  to  this  process  of  likely  false 
alarms,  each  contribution  to  the  objective  function  can  be 
weighted  in  accordance  with  the  table.  Choosing  actual 
numerical  weights  is  a  largely  empirical  process  based  on 
cut-and-try  procedures;  values  can  be  specified  and  adjusted 
by  the  user. 

A  related  consideration  is  the  difference  in  relative 
accuracies  of  x  and  y  components  of  detection  location  data. 
Recall  that  x  is  determined  by  position  along  a 
one-dimensional  CCD  array;  whereas  the  y  coordinate  is 
deduced  from  position  in  the  rotational  cycle  of  a  scanning 
mirror.  The  accuracies  of  these  two  physically  distinct 
processes  are  different.  Because  x  is  the  more  accurate, 
this  coordinate  should  be  accorded  more  effect  in  the 
optimization  process.  This  goal  is  accomplished  by 
asymmetrical  weighting  of  the  x  and  y  contributions  to  the 
objective  function.  Similarly,  the  dimensioning  of  the 
acceptance  cell  at  each  detection  point  must  accommodate  the 
greater  uncertainty  in  the  y-direction.  Accordingly,  the  box 
is  rectangular,  with  its  y-dimension  twice  that  of  x.  This 
asymmetry  can  be  seen  in  the  acceptance  cells  which  are  shown 
in  Figure  8.1. 


8 . 2 . 6 . 3  Diminishing  Size  of  Acceptance  Cells 

If,  among  n  identified  star-detection  associations,  nl  are 
caused  by  false  alarms,  it  can  be  expected  that,  whenever 
nx/n  is  sufficiently  small,  the  result  of  the  optimization 
process  will  be  close  star-detection  spacings  for  each  of  the 
individual  fn-n^  associations  corresponding  to  the  real 
detections;  but  the  spacings  for  the  nj  cases  of  false  alarms 
will  tend  to  be  substantially  greater.  It  is  desirable  to 
window  out  the  false  alarms,  thereby  enabling  a  more  accurate 
determination  of  attitude.  A  procedure  to  accomplish  this 
employs  multiple  stages  of  optimization.  At  successive 


stages  the  acceptance  window  is  narrowed  in  both  the  x  and  y 
directions.  The  concept  is  that  valid  star-detection  pairs 
will  be  relatively  closely  spaced  after  the  first 
optimization  and  will  fit  within  the  new,  smaller  cell. 

False  detection  pairs  will  be  more  widely  spaced  and  will  be 
rejected  by  the  smaller  cell.  Optimization  can  now  be 
repeated,  uncorrupted  by  the  false  data.  The  procedure  must 
be  truncated  whenever  (n-ni)  is  too  small  to  support  a 
meaningful  determination  of  attitude. 

8. 2. 6. 4  Multiple  Stars  within  an  Acceptance  Cell 

In  practice  more  than  one  star  may  fall  within  an  acceptance 
cell  centered  on  a  given  detection.  The  choice  among  the 
candidate  stars  is  made  on  the  basis  of  which  of  them 
corresponds  to  the  greatest  X25  optical  intensity. 
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